Mathematical Programming Computation
https://doi.org/10.1007/s12532-020-00194-3
FULL LENGTH PAPERMIPLIB 2017: data-driven compilationof the 6th mixed-integer programming library
Ambros Gleixner, Gregor Hendel, Gerald Gamrath, et al. [full author details at
the end of the article]
Received: 25 July 2019 / Accepted: 10 September 2020
© The Author(s) 2021

Abstract
We report on the selection process leading to the sixth version of the Mixed Integer
Programming Library, MIPLIB 2017. Selected from an initial pool of 5721 instances,
the new MIPLIB 2017 collection consists of 1065 instances. A subset of 240 instances
was specially selected for benchmarking solver performance. For the first time, thesesets were compiled using a data-driven selection process supported by the solution ofa sequence of mixed integer optimization problems, which encode requirements ondiversity and balancedness with respect to instance features and performance data.
Keywords Mixed integer linear optimization · MIP · Benchmarking · Selection
methodology · Instance library
Mathematics Subject Classification 90C06 · 90C09 · 90C10 · 90C11

1 Introduction
Computational mixed integer (linear) optimization is an important sub-field of mathematical optimization. Hundreds of papers on the subject are published each year anda multitude of companies provide tools for modeling and solution of mixed integer
optimization problems (MIPs) based on state-of-the-art techniques. Measuring performance on benchmark test instances has lain at the heart of computational research
since the early days of mathematical optimization. Hoffman et al. [28] first reported
on a computational experiment comparing implementations of three algorithms forlinear optimization back in 1953. Their observation that “[many] conjectures about
the relative merits of the three methods by various criteria could only be verified byactual trial” seems to hold to an even greater extent today. The variety of different

B Ambros Gleixner
gleixner@zib.de

Extended author information available on the last page of the article

123

A. Gleixner, G. Hendel, G. Gamrath et al.

techniques and their complex interactions regularly calls for empirical evaluation andmotivates the collection and curation of relevant test instances.
Brought into existence in 1992 by Bixby, Boyd, and Indovina [8], the goal of the
MIPLIB project has been to provide the research community with a curated set ofchallenging, real-world instances from academic and industrial applications that aresuitable for testing new algorithms and quantifying performance. It has previously
been updated four times [3,9,32,39] in order to reflect the increasing diversity andcomplexity of the MIPs arising in practice and the improved performance of available
MIP solvers. In this article, we describe its sixth version, MIPLIB 2017, together witha new selection methodology developed during this process.
The exceptional algorithmic progress in solving real-world MIP instances overthe last decades is recorded in various articles [7,32,33,35]. Specifically, this can beobserved by examining results on the previous version, MIPLIB 2010, both in termsof solvability and speed. By the end of 2018, the number of unsolved instances was
reduced by nearly half. Of the 134 instances for which no solution with provable
optimality guarantee was initially known, only 70 remain open. Comparable progressin the overall speed of solvers can be observed in the results of benchmark testingwith different versions of available solvers. Since its release in April 2011, the subset
of instances of MIPLIB 2010 that form the so-called “benchmark set”, consisting of
87 problem instances, has been the accepted standard for evaluating solvers. Using
this benchmark set, Hans Mittelmann has been evaluating a number of MIP solvers,
including CPLEX [13], Gurobi [25], and Xpress [51]. When MIPLIB 2010 was
released, the version numbers of these three commercial solvers were CPLEX 12.2.0.2,
Gurobi 4.5.1, and Xpress 7.2. Aggregating the benchmark results of these threesolvers at that time, we can construct results corresponding to a so-called “virtual
best” solver and a so-called “virtual worst” solver. These are hypothetical solvers that,for each instance, produce run times that are equal to the best and the worst of the three,respectively. Doing this analysis yields shifted geometric mean runtimes of 36.3 and113.0 s for the virtual best and virtual worst solver, respectively.1 In December 2018,the solver versions were CPLEX 12.8.0, Gurobi 8.1.0, and Xpress 8.5.1. On the same
hardware (with a newer operating system) the shifted geometric means of the runtimes
had decreased to 13.5 s for the virtual best, and 31.3 s for the virtual worst solver.
This corresponds to speed-up factor of 2.70 and 3.62, respectively, which amounts toroughly 16% per year, just from improvements in the algorithms.
It was because of this development that the MIPLIB 2010 benchmark set was nolonger considered to sufficiently reflect the frontier of new challenges in the field andthe process of constructing a new MIPLIB began. In November 2016, a public callfor contributions was launched and a group of 21 interested researchers, including
representatives of the development teams of nine MIP solvers formed a committee inorder to steer the process of compiling an updated library.2 As with MIPLIB 2010,
1 The computations used 12 parallel threads. The corresponding log files can be found at [40]. The means
were computed with a shift of 1 second (see Achterberg [1]).
2 The members of the MIPLIB 2017 committee were Tobias Achterberg, Michael Bastubbe, Timo Berthold,Philipp Christophel, Mary Felenon, Koichi Fujii, Gerald Gamrath, Ambros Gleixner, Gregor Hendel, Kati
Jarck, Thorsten Koch, Jeff Linderoth, Marco Lübbecke, Hans Mittelmann, Derya Ozyurt, Imre Pólik, Ted
Ralphs, Domenico Salvagnin, Yuji Shinano, Franz Wesselmann, and Michael Winkler.

123

MIPLIB 2017: data-driven compilation...

the overall goal was the compilation of two sets of instances. The MIPLIB 2017
benchmark set was to be suitable, to the extent possible, for performing a meaningful
and fair comparison of the average performance of MIP solvers (and different versionsof the same solver) across a wide range of instances with different properties, ina reasonable amount of computing time. The larger MIPLIB 2017 collection was
to provide additional instances for a broader coverage without restrictions on the
total runtime of the test set, including unsolved instances (as a challenge for future
research) and instances not suitable for benchmarking due to problematic numericalproperties, special constraint types (such as indicators), or exceptionally large memory
requirements.
It should be emphasized that the benchmark set we have scientifically constructed
is designed for the purpose of comparing the overall performance of general purpose
solvers on a wide-ranging set of instances. Average performance on this set is not asuitable criterion to decide which MIP solver to use in a particular application scenario.
For such decisions, it is important to consider what specific class(es) of instances are
relevant, as well as what criteria beyond the raw speed and the ability to solve a wide
range of problems are of interest. This is also underlined by the fact that each of theeight solvers that were used to collect performance data (see Sect. 4.6) proved to be
the fastest solver on at least one instance.
Compiling a representative and meaningful instance library is a nontrivial endeavor.
Compared to previous editions of MIPLIB, the increased number of submissions, thegoals of compiling a significantly larger collection of instances and including a larger
number of representatives of solvers posed new challenges to the selection process. In
addition, MIPLIB 2017 is the first edition to provide supplementary data regarding the
instances, such as the matrix structure and decomposability, as well as the underlying
models from which the instances originated, where available. In order to produce awell-balanced library in a fair and transparent manner, we designed a new, heavily
data-driven process. The steps applied between the initial submissions and the final
MIPLIB 2017 are outlined in Fig. 1. Driven by a diverse set of instance features, ourmethodology used multiple clusterings to populate a MIP model that was then solved
to generate suitable candidates for the final library to be presented to the MIPLIB
committee.
We consider this process of selecting instances from a large pool of submissions tobe the main new feature of MIPLIB 2017. By contrast, the instances constituting previous versions of MIPLIB were manually selected by the members of the committee,
depending heavily on their expertise in benchmarking to avoid common pitfalls likeoverrepresentation of certain problem classes. As one byproduct of this data-driven
approach, we are now able to identify similar instances, which leads to sometimes
surprising insights into connections between different, seemingly unrelated instancesin the library. In addition to the raw feature data, we provide, for each instance, the five
most similar instances in the collection on a newly designed web page (see Sect. 6.3).
Despite the heavy use of quantitative methodology, the overall process inherently
requires many manual decisions and heuristic choices. Some examples are the choiceof features, the number of clusters to use when clustering instances according to thisfeature data, the definition of which instances to consider as computationally easyor hard, and our formalizations of diversity and balancedness with respect to the123

A. Gleixner, G. Hendel, G. Gamrath et al.

Fig. 1 Outline of the steps involved in the selection of the MIPLIB 2017 collection and benchmark set.
Number of instances remaining are given in parenthesesfeature data. All of these decisions have a high impact on the final result. Hence, in theremainder of this article we try to describe the collection and selection of MIPLIB 2017at a sufficiently high level of technical detail in order to make this process transparent
to the reader.
The article is organized as follows. Section 2 addresses related work and gives abrief historic overview on standards for computational experiments in mathematicaloptimization. Section 3 describes the efforts to collect instances and meta data on theirbackground and solvability. Section 4 details the collection of feature data for each
instance that forms the basis for quantifying similarity of instances and balancednessof a given selection. In Sect. 5, we describe how this data was used as input in order to123

MIPLIB 2017: data-driven compilation...

compute good candidates for the library by solving a MIP model. Section 6 summarizes
the final result. We conclude with our final remarks and acknowledgements in Sect. 7.

2 Related work
With the evolution of computational research, standards and guidelines for conducting
computational experiments were proposed and updated. Next to performance measures, software availability and machine specifications, these guidelines also discuss
the choice of test problems and the availability of standard test sets.
In 1978, Crowder et al. [14] presented an early list of best practices, emphasizing
the necessity of reproducibility of computational experiments. In this course, they also
discussed the value of test problems being documented and the importance of using
real-world problems. At the same time, they pointed out that sharing test problems is
expensive and time-consuming—an issue that we have fortunately overcome.
The work of Jackson et al. [31] was motivated by a controversy around the firstpublished computational results for interior point LP solvers. Besides others, the adhoc committee stated: “We recommend that standard test problem sets and standard test
problem generators be used wherever possible.” Further, they encouraged researchersthat “whenever new problems are used in a computational evaluation, these problems
be …submit[ed] to a standard collection like netlib.”
In his seminal paper from 1993, John Hooker characterized an empirical scienceof algorithms [29]. He discussed the importance of identifying typical representativesof a problem class to conduct a study on, at the same time mentioning that any choice
is open to the criticism of being unrepresentative. There are several resolutions forthis issue; a well-known, long-standing and publicly available standard test set iscertainly one of them. Another resolution that Hooker pointed out is to make “the issue
of problem choice …one of experimental design”. This means making the question
how performance depends on test set characteristics part of the experiment itself.
The approach taken to set up MIPLIB 2017—a data-driven instance selection thatparameterizes the creation of a test set—can be seen as an extension of this idea.
There are various publications that formalize the problem of compiling a test set.
McGeoch [37,38] developed an abstract model of algorithms and the paradigm ofsimulation research, partially in response to Hooker’s paper.
A complementary line of research that is not touched on in this paper is the creation of new test instances to fill “gaps” in a test set, by learning from the instanceparameter space, see [49]. Smith-Miles et al. [48] use such an approach to work out
the strengths and weaknesses of optimization algorithms on an enlarged instance set
that they extrapolated from a standard test set. They detect so-called pockets wherealgorithm performance significantly differs from average performance. This take on
instance diversity complements our approach of trying to achieve a best possible coverage of the feature space subject to a fixed set of candidate instances.
Work on standard test sets for benchmarking naturally connects to work on algorithm selection [43]. The work of Bischl et al. [6] brings both fields together by
publishing a benchmark library for algorithm selection. They introduce a data format

123

A. Gleixner, G. Hendel, G. Gamrath et al.

to formally define features, scenarios, working limits and performance metrics. Since
2016, this library contains a MIP section, based on MIPLIB 2010.
Finally, related work of course includes various other libraries for mathematicaloptimization problems, including, but not limited to MINLPLib [11], QPLIB [16],Netlib [10], the COR@L collection [34], and OR-Library [4]. The latter is one of theoldest online collections of optimization problems, existing for 30 years now and still
being regularly updated.

3 The collection process
The first step in the creation of MIPLIB 2017 was to collect a large set of candidate
instances. This section elaborates on the submission process and the origin of thesubmitted instances, as well as preparatory steps for the subsequent selection process,such as instance cleanup and grouping.
3.1 The submission phase
The collection process started with a call for contributions in November 2016, which
was advertised at many conferences, on numerous mailing lists, and in direct communication with partners from industry. Submissions were accepted until August 1, 2017.
Overall, we received 128 submissions (each containing multiple instances) from 66
submitters with 54 different affiliations, 38 of them being academic and 16 industrial.
Note that the affiliation of the submitter alone does not imply anything about the sourceor type of the model. Several of the submitters with academic affiliation submitted
instances from industry projects and one submitter from industry submitted instancesmodeling a combinatorial game.
Each submission was automatically committed as a new subdirectory to the public git repository https://git.zib.de/miplib2017/submissions. In the repository, thesesubdirectories are structured as follows. In the root of the subdirectory there are two
files: a bibtex file with references to papers related to the instances, if provided by
the submitter, and a meta file containing information about the submitter, creator, andowner of the submitted instances, as well as a description of the instances and licensing
information. For the first time, all contributions were required to be submitted under alicense in order to explicitly grant the right of redistribution, among others. The default
license was chosen to be the Creative Commons CC BY-SA 4.0.3 license, but it waspossible to specify a different license if desired. In addition to these two files, there
are up to three subdirectories. The instances subdirectory contains the instances
themselves (.lp or .mps format, possibly compressed). The models subdirectory contains associated model files, which contributors were encouraged to include in order toprovide researchers with richer information on the structure of the instances. Finally,
additional information is provided in the misc subdirectory, ranging from extended
descriptions to log files and MIP start files.
3 https://creativecommons.org/licenses/by-sa/4.0/.

123

MIPLIB 2017: data-driven compilation...

In total, 3670 instances were newly submitted to MIPLIB 2017. Table 1 lists all
submitters and the number of instances they submitted. We arrived at a total of 5721
instances by adding 2051 instances that were submitted for inclusion in MIPLIB 2010,
keeping their original submission information intact. Those 2051 comprised most of
the submissions to MIPLIB 2010 except for a few duplicates already present in other
MIPLIB 2017 submissions.
Table 2 reports the origin of the submitted instances. It shows that there are two
blocks of new submissions with a high number of instances: the instances submitted
to the NEOS server and MIP models of instances that were part of the MiniZinc
Challenges4 from 2012 to 2016. Even excluding those instances, the remaining newcontributions comprise about as many instances as were submitted to MIPLIB 2010,
half of which were collected from publicly available sources back then. Note also thatnearly 50 instances originate from submissions to the user support of two commercial
solvers.
3.2 Submissions from the NEOS Server
Before moving on, we briefly discuss how the instances originating from the NEOS
Server were collected, since this was a separate procedure unto itself. The NEOS
Server is a free internet-based service for solving numerical optimization problems
hosted by the Wisconsin Institute for Discovery at the University of Wisconsin in
Madison, with remote solving services provided by various sites, such as ArizonaState University.
In the calendar years 2014–2016, more than 1 million MIP instances were solved
using NEOS. A subset of these instances was collected and submitted for inclusion into
MIPLIB 2017. To avoid misunderstandings, we note that the reduction methodology
described in the following happened prior to the MIPLIB 2017 effort and is not partof the selection procedure conducted for MIPLIB 2017. However, because the NEOS
instances account for a large number of newly submitted instances, we outline theprocess briefly for the sake of completeness and transparency.
To begin, the 14,571 MIP instances that were modeled through the AMPL [15]language and whose solution required more than five minutes of wall clock time
to solve were collected. In discussions with the MIPLIB committee, it was decided
that a representative, yet diverse, sample of around 700-800 instances would be an
appropriate number of submissions from NEOS. The strategy to select the instances
was based on clustering the instances with similar traits and then selecting instancesfrom the clusters.
All 14,571 of these instances were re-run using CPLEX 12.6 to collect the optimal
solution (if possible within 12 h) and the solution to the root LP relaxation. For each
instance, the following properties were collected:– the number of variables, constraints, and nonzeros,– the percentage of binary variables and of general integer variables,– the best solution value found by CPLEX, the root gap, and the percentage of binary
or integer variables with fractional value in the root LP solution.
4 https://www.minizinc.org/challenge.html.

123

A. Gleixner, G. Hendel, G. Gamrath et al.
Table 1 Submitters of new instances, and number of instances from their submission(s) after cleanup
Submitter#

Submitter#

Andrea Arias

360

Felix Cebulla

Dimitri Papageorgiou

300

Irv Lustig

3

Pierre Le Bodic

121

Jesus Rodriguez

3

Qie He120

Jeff Linderoth (NEOS server)

Simon Felix

101

Jeff Linderoth (other)

3

Gleb Belov (MiniZinc challenge)

988

Jonathan Eckstein

3

Gleb Belov (other)100

Siwei Sun

3

Simon Bowly

92

Felix J L Willamowski

2

Gerald Gamrath76

Janos Hoener

2

Hans Mittelmann (NEOS server)

75

Joshua Friedman

2

Haroldo Gambini Santos

71

Utz-Uwe Haus

2

Shunji Umetani

60

Andreas Baermann

1

Matias Soerensen

42

Balabhaskar Balasundaram

1

Jordi Castro

34

Christian Liebchen

1

Toni Sorrell

34

Christopher Daniel Richards

1

Stephan Beyer

32

Dan Neiman

1

Manuel Iori

30

Daniel Bienstock

1

Michael Winkler

28

Daniel Rehfeldt

1

Cezar Augusto Nascimento e Silva

26

Gavin Goodall

1

3

710

Pelin Damci-Kurt

26

Gerald Lach1

Domenico Salvagnin

22

Hsiang-Yun Wu

1

Michael Bastubbe

19

Jesse Liu

1

George Fonseca

18

Juan Javier Dominguez Moreno

11Sascha Kurz

17

Koichi Fujii

Antonio Frangioni

15

Mark Husted

1

Daniel Heinlein

12

Paula Carroll

1

Marc Pfetsch

11

Sujayandra Vaddagiri

1Tamas Terlaky

11

Timo Berthold

1

Berk Ustun

10

Philipp Leise

9

Salim Haddadi

9

Dan Hiroshige, Koichi Fujii

8

Christopher Hojny

7

Laurent Sorber

6

Yoshihiro Kanno

6

Andrew Stamps

5

Alexandra M Newman

4

Austin Buchanan

4

Rob Pratt4

Sean MacDermant

4123

MIPLIB 2017: data-driven compilation...
Table 2 Origin of the submitted
instances. Each instance iscounted only once, in the first
applicable row

MIPLIB 3

65

MIPLIB 2003

28

MIPLIB 2010

335

MiniZinc Challenge 2012–2016

988

New submissions (NEOS server)

785

New submissions (other)

1897

MIPLIB 2010 submissions

1623

Instances whose number of variables or constraints was more than 3 standard deviations away from the mean were excluded as outliers. In an attempt to eliminate
duplicate instances, only one instance in a group that had the same number of constraints, same number of variables, and same optimal solution value, as well as the same
associated email address, was retained. After this removal of outliers and suspected
duplicates, 6531 instances remained.
Each of these instances was assigned to a multi-labeled category using three independent calls to a k-means clustering algorithm [26] as follows. First, the instanceswere divided into four clusters according to properties associated with problem size:
the number of constraints, number of variables, and number of nonzeros. Second, theinstances were again divided into four clusters, this time according to the percentage
of binary and general integer variables in the instance. Finally, the instances weredivided into four more clusters based on two properties of the root LP relaxation: root
gap and the percentage of integer variables whose root LP relaxation value was not
integer-valued.
From these clusterings, each of the 6531 instances was given a label in {0, 1, 2, 3}3 ,all in all giving 43 = 64 (possibly empty) clusters. Let Si be the number of instancesin the i-th cluster. The final selections were then made as follows:– For each cluster with size Si ≤ 2, all instances of the cluster were selected;– for each cluster with size 3 ≤ Si ≤ 10, |Si |/2 instances were selected at random;– for each cluster with size Si > 10, |Si |/10 instances were selected at random.
This strategy ensured a diverse sample, and resulted in a total of 710 instances to beconsidered for inclusion in MIPLIB 2017.
Separately, Hans Mittelmann also contributed 75 instances submitted to one of the
ASU-hosted servers through the NEOS platform. These servers process requests for
the SCIP, FEASPUMP, PROXY, and Qsopt_ex solvers. Using scripts, these user
submissions were screened for promising instances with an emphasis on those that
contain model information.
In combination, this resulted in a total of 785 new instances arising from NEOS. In
addition, the instance pool for MIPLIB 2017 contains 391 NEOS instances that havebeen either included in previous versions of MIPLIB or were considered during past
selection processes.

123

A. Gleixner, G. Hendel, G. Gamrath et al.

3.3 Instance cleanup
Following the close of submissions, we performed several cleanup steps, as follows.
Instance renaming While we tried to preserve original file names where possible, some submissions required renaming to create uniquely identifiable names. For
instance, some submissions contained multiple instances with the same file name thatoriginated from the same base model populated with different data. We also tried tomake the names as short and expressive as possible, containing only lower-case characters, numbers and hyphens. After re-naming, the longest instance name was reduced
from 151 to 31 characters. Note that submissions for MIPLIB 2010 were renamed tobe consistent with the naming scheme of the final MIPLIB 2010 instances, but nosubmissions already present in MIPLIB 2010 were renamed.
Making NEOS instances identifiable NEOS instances have traditionally been named
neos-xyz, where xyz is a number with up to 7 digits, representing the unique NEOS
instance ID of the submitted instance. In order to allow easier identification of theinstances in papers and personal communications, we compiled a list of river names
from all over the world and appended a river name to the original name, resulting
in names of the form neos-xyz-river. As an example, NEOS instance 3075395 has
been renamed neos-3075395-nile such that it can be colloquially called the “nileinstance”. We excluded river names such as “Amazon” to avoid ambiguous renaming.
Note that we applied this renaming procedure only to the 785 NEOS instances newly
submitted for MIPLIB 2017 (see Sect. 3.2), leaving all previously available NEOS
instances under their old name to avoid confusion.
Format conversion and cleanup All instances in the MIPLIB 2017 collection areprovided in MPS format [30,42]. Therefore, instances submitted in LP format wereread into CPLEX and written out in MPS format. Given that different MIP solvers
support different extensions of the MPS format, we ensured that all solvers could
read all instances by restricting the format to a basic version. Maximization problems
were turned into minimization ones by negating objective coefficients; lazy constraint
markers were removed so that the constraints are treated as standard constraints; andcoefficients stated as hexadecimal numbers were converted to their decimal equivalent.
Additionally, we added a NAME section or changed the existing one to the updated
instance name. For a small number of instances, it was necessary to change the MPS
files because they contained ambiguous definitions or outright incorrect formatting.
In those cases, we performed the minimally necessary changes to make the instancefile adhere to the basic MPS standard.
3.4 Model groups
In most cases, the instances in a single submission are closely related in the sense thatthey originate from the same or a very similar MIP model. Some submissions, however,contain many unrelated instances. Therefore, we introduced model groups to keep track

123

MIPLIB 2017: data-driven compilation...
Table 3 Model groups and counts for the different instance sources. Each group is counted only in the first
applicable row
TypeGroupsGroup size ∈{1}{2, . . . , 5}

{6, . . . , 10}{11, . . . , 360}

MiniZinc

80

0

16

38

26

NEOS

110

8

24

34

44

MIPLIB 2017 submissions

130

81

15

10

24

MIPLIB 2010 submissions

241

172

30

16

23

of this form of meta information that may not be directly inferable from the submission
ID or the numerical instance features described in Sect. 4. A model group represents a
set of instances that is based on the same model or a very similar formulation but withdifferent data. This grouping allowed us to avoid overrepresentation of a particularapplication or model class in the final library by limiting the number of instances withknown similar model background during the selection process.
Each instance was assigned to one model group as follows. Initially, a submissionof homogeneous instances was assigned to its own model group. If a submission
contained multiple sets of instances, each implementing a different model for the sameproblem and data set, an individual group was created for each of the different modeltypes. This procedure was applied to both the submissions to MIPLIB 2017 and thesubmissions to MIPLIB 2010. Publicly available instances with known application
were grouped by hand by the authors. Examples for such cases are the MiniZinc
instances submitted to MIPLIB 2017 and the instances from older MIPLIB versions
and public instance sets submitted to MIPLIB 2010.
Instances from the NEOS server, however, are anonymous and lack meta data fromthe submitters that could be used to infer model groups. Nevertheless, users often
submit multiple, similar instances. Hence, we used feature data described in the next
section in order to infer synthetic model groups in an automated way. In order to group
the NEOS instances, a k-means clustering was computed with respect to the entire
instance feature space (see Sect. 4). The parameter k = 110 was chosen manuallyto achieve a clustering with very similar instances in each NEOS model group. This
clustering was applied to all 1176 NEOS instances both from new submissions andpreviously available sources.
Table 3 summarizes the number of resulting model groups and the corresponding
group sizes for the different sources of instances. These numbers are given with respect
to the 5666 instances in I sub (see Fig. 1). The largest model group cmflsp comprises
360 instances of a capacitated multi-family lot-sizing problem.

4 Feature computation
The ultimate goal of the selection process was to select a representative sample of
all available instances with respect to problem structure and computational difficulty,
while avoiding overrepresentation of similar models. In the spirit of data-driven deci-

123

A. Gleixner, G. Hendel, G. Gamrath et al.

sion making and in light of related work in the fields of algorithm selection and machinelearning, we based this process both on performance data and on an extensive set of
instance features. The first step in the selection process was simply to determine thefeatures of interest. In the terminology of [48], this means that we defined a feature
space of measurable characteristics and computed a feature vector associated to eachelement of the problem space of candidate instances.
Although this may seem straightforward, it is important to note that the feature
vector corresponding to an instance can be affected by seemingly irrelevant propertiesof its representation in MPS format. For instance, some of the raw MPS instances
contained modeling language artifacts or artifical redundancies. For this reason, theinstance features were computed only after applying some straightforward simplification steps, which we refer to as trivial presolving. We first describe this presolving
process before describing what features of the instances were used and how theirvalues were determined for each instance.
4.1 Trivial presolving
As is traditional for MIPLIB, the submitted instances are being made available in their
original, unmodified form (with the exception of minor corrections to MPS formattingthat were necessary in some cases). This means that the distributed instances were notpresolved or simplified in any way for the purposes of distribution. For the purposesof extracting features of the instances, however, so-called “trivial” presolving wasapplied, as described below. It may seem strange that the version of each instancemade publicly available in the final collection may actually be slightly different thanthe version considered during the selection process, but there are good reasons for thisapproach that we elaborate on next.
The justification for distributing the instances in their original submitted form is
simply that this allows the most complete and realistic testing of the ability of eachsolver to deal with real-world instances, including all of the idiosyncratic artifacts that
may arise in the modeling process. In particular, algorithms for presolving instances areactively developed and have a high impact on the performance of a solver (see [2,20]).
They not only strengthen the original formulation, but also simplify and remove unnecessary artifacts. These procedures are computational in nature and their efficiency andeffectiveness also needs testing. In some cases, there may be choices made in thepresolving that can only be made with foreknowledge of the properties of the solution
algorithm itself. For all these reasons, it was considered highly desirable in promoting test conditions that are as reflective of real-world conditions as possible to avoid
modifying the distributed versions of the instances.
On the other hand, because MIP solvers do universally apply certain well-known
simplification procedures to an instance before the branch-and-bound search, theunmodified descriptive data of the original instance may not properly reflect the “true”
features of that instance for the purposes of clustering instances according to similarity,
as we did during the selection process. The features considered that may be affected bypresolving include not only obvious properties, such as instance size, but less obvious
ones, such as the type of constraints and variables present in the model. A model may,123

MIPLIB 2017: data-driven compilation...

for example, have all variables integer with the exception of one continuous variable
whose value is fixed to 1 and whose purpose is to model an objective offset. It would
be unreasonable to consider such an instance to be an instance with both integer andcontinuous variables. At the other extreme, we may have an instance in which allbinary and integer variables are implicitly fixed, leaving a purely continuous problem
after presolving.
While it seems necessary to do some presolving before computing instance features,the full presolving done by solvers is itself a difficult computational balancing act andeach solver does it differently. Too much presolving before feature computation would
result in a presolved instance with features no more representative of the “true” ones
than the completely unpresolved instance. As a compromise, all instance features
introduced in Sects. 4.3–4.5 were collected after applying a reduced set of the mostobvious presolving techniques to the instance, but no more sophisticated techniques.
For this trivial presolving, we used SCIP 5.0, but disabled most presolving techniques, applying only simple ones, such as the removal of redundant constraints andfixed variables, activity-based bound tightening, and coefficient tightening. In contrast
to standard SCIP presolving, which stops if the problem size could be reduced by onlya small percentage during the last presolving round, we applied the simple presolving
steps until a fixed point was reached. The complete set of SCIP parameters used to
do the presolving is provided on the MIPLIB web page (see Sect. 6.3) as part of thefeature extractor download.
For 55 of the 5721 submitted instances, trivial presolving turns the instance into apure LP or is even able to solve the instance by fixing all variables. These instanceswere not considered for inclusion and also serve to emphasize the importance of
this preprocessing step. Overall, trivial presolving reduced the number of variables on3782 instances (66% of the submission pool), sometimes by as much as 93% (instancea2864-99blp). For 445 instances (8%), more than 50% of the variables were fixed.
On average, trivial presolving reduced the number of variables by 15%.
4.2 Canonical form
Because the feature computation can be affected not only by presolving but also by theexact form in which the instance is represented (equality constraints versus inequalitiesetc.), we transformed all presolved instances into the following canonical form, whichis slightly more general than the usual one, prior to feature computation.
Definition 1 A mixed integer optimization problem P with input– m, n, n b , n i , n c ∈ N, n = n b + n i + n c ,– coefficient matrix A ∈ Qm×n ,– left-hand and right-hand side vectors  A , u A ∈ Qm±∞ ,– lower and upper bound vectors x , u x ∈ Qn±∞ , and– objective coefficient vector c ∈ Qn
is defined to be an optimization problem of the form
min cx :  A ≤ Ax ≤ u A , x ≤ x ≤ u x , x ∈ {0, 1}n b × Zn i × Qn c .

123

A. Gleixner, G. Hendel, G. Gamrath et al.

It is important to note that transformation to the canonical form of Definition 1 is
not uniquely determined for each instance. There remain certain degrees of freedom toformulate equivalent instances by scaling continuous columns, scaling the objectivefunction c x, or scaling a constraint iA ≤ ai x ≤ u iA . This may cause problems
with the computation of some features. For example, some features of interest involve
comparison of row coefficients, but this comparison is difficult if rows of the constraint
matrix may have coefficients differing by several orders of magnitude. We addressthese issues by normalizing the objective coefficients c and every constraint iA ≤ai x ≤ u iA by their maximum absolute coefficient c∞ and ai ∞ , respectively, sothat all objective and matrix coefficients lie in the interval [−1, 1] before computing
the feature matrix F (the downloadable instances are not altered).
It is based on this final presolved canonical representation that we define the Q =105 features we consider. This results in a feature matrix F ∈ R N ×Q , where N is thetotal number of instances submitted. Table 4 lists these features, which were divided
into K = 11 feature groups (also listed in the table). Feature groups were used forthe selection process, during which instance clusters were computed for each featuregroup individually. Every feature group was chosen to represent a particular aspect of
an instance in the form specified by Definition 1. The computation of features in mostof the groups only requires information that can be extracted directly from the inputof the (trivially presolved) problem. Two exceptions are the constraint classificationand decomposition groups, which need to identify structures in the model. These aredescribed in Sects. 4.4 and 4.5 .
4.3 Instance features
Here, we describe the first nine feature groups in Table 4. We use the shorthand vector
statistics to refer to five values summarizing the entries of a vector v ∈ Rd±∞ . Let
d = |{ j : |v j | < ∞}| be the number of finite entries of v, which can be smaller thand in the case of, e.g., bound vectors, and let v be the restriction of v to its finite entries.
We assume without loss of generality that v is sorted, v1 ≤ v2 ≤ · · · ≤ vd . The fivevalues are– min : v → v1 ,– max : v → vd ,– mean : v → d1

d


vj,



– median : v → v d +1 + v d +1 /2, and
 2 
2

d 
2

v j − mean(v ) .
– std : v → d1
j=1

j=1

Note that infinite entries can only occur for the variable bound vectors x and u x andthe left- and right-hand side vectors  A , u A . For a vector v that contains only infinite
entries, i.e., for which d = 0, the above vector summaries are not well-defined. If
d = 0, the corresponding statistics were set to 0 in the data. Note that even if the123

MIPLIB 2017: data-driven compilation...
Table 4 Description of instance features used. Set notation is abbreviated, e.g., {A  = 0} denotes {(i, j) ∈{1, . . . , m} × {1, . . . , n} : ai, j  = 0}Group

Features

Description

Scaling

Size

3Size m,n of matrix, nonzero entries|{A  = 0}|log10 (x)2

Variable
types

3

Proportion of binary, integer, andnncontinuous variables nb , ni , nnc

Objectivenonzerodensity5

Nonzero density of objective function|{c =0}|
both total and by variable typen(bin., int., cont.), 0–1 indicator for
feasibility problems without objective

Objective
coefficients

6

Vector statistics and dynamism of c

Variablebounds12

Matrixnonzeros

6

Vector statistics of nonzero entries|{ai  = 0}| by row in A, nonzeros percolumn |{ An=0}|log10 (x) for
nonzeros percolumnMatrixcoefficients

19

Vector statistics of the fourm-dimensional vectors describing themin, mean, max, and std of the nonzero
coefficients {ai  = 0} in each row

Every ainormalizedby ai ∞

Row
dynamism

5

Vector statistics of row dynamism

log10 (x)

Sides

19

Vector statistics of left- and right-handsides  A , u A and concatenated(| A |||u A |), nonzero and finite
densities of  A , u A

Every ainormalizedby ai ∞Constraintclassification

17

Proportion of classes of special linear
constraints: singleton, precedence,
knapsack, mixed binary (see Sect. 4.4)

Decomposition

10

Features describing decomposition D
found by GCG with maximum
area score (see Sect. 4.5):
areascore(D), k, vector statistics
 r
|Dkr ||D1 |(except std) of
m ,..., m c|Dkc |
|D1 |
and
,
.
.
.
,
. Not available
n
n

x

x

Finite densities |{| |<∞}|, |{|u |<∞}|nnof bounds, vector statistics of upperx
bounds u and bound ranges u x − x .

ai ∞min j {|ai j | =0}

c normalizedby c∞Vectorstatistics
scaled bysiglog(x)for all instances123

A. Gleixner, G. Hendel, G. Gamrath et al.

original formulation has infinite bounds on variables, trivial presolving may often infer
finite bounds for those variables.
The dynamism of a vector with finite entries is the ratio of the largest and smallest
absolute entries, i.e., v∞ /min{|v j | : v j = 0}. The dynamism is always at least 1. If
the dynamism of any single constraint exceeds 106 , this is an indication of a numerically difficult formulation. Note that the dynamism is invariant to the normalization
procedure. Combining the dynamism of each constraint yields an m-dimensional vector, which can be summarized using vector statistics.
The feature group Matrix coefficients summarizes the nonzero coefficientsof the matrix A as follows. First, each row ai , 1 ≤ i ≤ m, of A is normalized byits largest absolute coefficient, such that all coefficients are in the range [−1, 1]. The
nonzero entries of ai are then summarized by four of the five vector statistics explained
above, namely the min, max, mean, and std. Going through all rows, we obtain fourm-dimensional vectors describing the min, max, mean, and std per row. Each of these
vectors is then summarized via vector statistics, which yields a total of 20 statistics
that summarize the cofficients of A. Examples are the mean minimum coefficient over
all m rows, or the standard deviation of all m maximum coefficients, etc. The featuregroup comprises 19 out of these 20 coefficient statistics, because the maximum over
all m maximum coefficients is equal to 1 for every instance in our data set.
For the feature group Sides, the m-dimensional left- and right-hand side vectors  A
and u A are summarized individually via vector statistics of all their finite elements.
Besides, we compute vector statistics for the finite elements of the concatenated 2mdimensional vector (| A |, |u A |) that combines the absolute left- and right-hand sidesof all rows. Note that the row normalization by the maximum absolute coefficient also
affects the row left- and right-hand sides.
For features such as the row or objective dynamism, which may differ by ordersof magnitude between instances, we used a logarithmic scaling. While logarithmic
scaling is fine for vectors with positive entries, it is not applicable to vectors withpotentially negative entries such as the variable lower and upper bound vectors. In
those cases, we apply a customized scaling
siglog : R → R, x → sig(x) log10 (|x| + 1)
to every entry of the corresponding column in the feature matrix F. The map siglog
preserves the sign of each entry.
The collection of the instance features was performed with a small C++ application
called the feature extractor, which extends SCIP by the necessary functionality needed
to report features after trivial presolving and optionally accepts a settings file to modify
the default presolving explained in Sect. 4.1. The feature extractor is a modified version
of a code used already by [22] and available for download on the MIPLIB 2017 web
page (see Sect. 6.3). 5
5 The actual computations reported in the following were carried out with five additional, redundant matrix
features. Only during the preparation of the manuscript, they were identified to be identical to other features.

123

MIPLIB 2017: data-driven compilation...
Table 5 Classification of linear constraints, sorted from most specific to most general. A constraint isalways assigned the first (topmost) type that applies. Notationally, a, b, c denote coefficients and righthand, possibly restricted to integer and/or nonnegative values, while x, y, and z are variables
Type

Linear constraints …

Empty

with no variables

Freewith no finite side

Singletonwith a single variable

Aggregationof the form ax + by = c

Precedenceof the form ax − ay ≤ b where x and y are binary/integer/continuous

Mixed Binaryof the form ax + by ≤ c, x ∈ {0, 1}

xi = 1, xi ∈ {0, 1} ∀i
of the form xi ≤ 1, xi ∈ {0, 1} ∀i
of the form xi ≥ 1, xi ∈ {0, 1} ∀i
of the form xi = k, xi ∈ {0, 1} ∀i, k ≥ 2
of the form xi ≤ b, xi ∈ {0, 1} ∀i, b ∈ N, b ≥ 2
of the form ai xi = b, xi ∈ {0, 1} ∀i, b ∈ N, b ≥ 2
of the form ai xi + ax ≤ a, x, xi ∈ {0, 1} ∀i, a, ai ∈ N∀i, a ≥ 2
of the form ai xi ≤ b, xi ∈ {0, 1} ∀i, b ∈ N, b ≥ 2
of the form ai xi ≤ b, xi ∈ Z ∀i, b ∈ Nof the form ai xi + b j z j {≤, =} c, xi ∈ {0, 1} ∀i, z j ∈ R ∀ jGeneral Linearwith no special structure

Variable Bound
Set Partitioning
Set Packing
Set Covering
Cardinality
Invariant Knapsack
Equation Knapsack
Binpacking
Knapsack
Integer Knapsackof the form

4.4 Constraint classification
Table 5 lists the constraint classification types used for the feature group Constraintclassification. A total of 17 types of linear constraints that often occur as a subsetof the constraints of MIP instances were identified. The table is sorted from most
specific to most general. If a constraint belongs to multiple types, the classification
always assigns the most specific, i.e., topmost, type that applies. Note that even empty,
free, and singleton constraints are listed. While these types are removed during trivial
presolving, they may well be present in the original formulation.
There are several types of constraints supported by the MPS format [30,42] that are
not strictly linear as required by Definition 1. A well-known extension are indicator
constraints, which are conditional, linear constraints that only need to be satisfied if acorresponding binary variable, the so-called indicator variable, is set to 1. It is possible
to linearize such a constraint by employing a sufficiently large coefficient M for theindicator variable, in which case the reformulation is called a big-M formulation.
In many practical applications, big-M formulations require a very large value of M,
which is why they often lead to numerically difficult models. Directly expressing
such constraints as indicator constraint allows the solver to handle them in a morealgorithmically advantageous way.
Indicator constraints were allowed into the MIPLIB 2017 collection, but (as we
describe later) were not allowed in the benchmark set. The only feature used that123

A. Gleixner, G. Hendel, G. Gamrath et al.

involves indicator constraints was their fraction among all constraints. This feature isalso part of the feature group Constraint classification in Table 4. Other features
regarding, e.g., the linear part of the indicator constraint, are not collected. In total,
437 of the submitted instances contain indicator constraints. Many of them appear ininstances of the MiniZinc submission, which were additionally submitted as big-M
formulations.
There are other special types of constraints allowed by MPS, such as special-orderedsets (SOSs), semicontinuous variables, and piecewise linear constraints that not all
solvers support. However, none of the instances submitted used such constraints.
4.5 Block-structured decomposition features
Dantzig-Wolfe reformulation and Lagrangian relaxation are decomposition techniquesthat can exploit the presence of logical groupings of constraints and variables in themodel. Classically, one identifies “complicating” or “linking” constraints that, whenremoved from the model, result in several subproblems that can be solved independently. Concretely, this occurs when the constraint matrix has so-called block angularstructure (defined below). One challenge in applying these techniques within a general purpose solver is identifying such structure. Standard input formats do not allow
information about these structures to be passed to the solver directly. For this reason,
several decomposition-based solvers have been developed that accept auxiliary input
files indicating this structure when it is present. These include GCG [21], DIP [17–19], and DECOMP (a decomposition-based solver that is part of SAS/OR and SAS
Optimization [44]). All three can exploit the identified structure by reformulating the
instance and applying a branch-and-price algorithm to solve it.
Although block structure may be “obvious” in the original model, this structure
is often lost when the model is populated with data and transformed into a format
such as MPS. When information about the block structure is not provided, it can still
be derived algorithmically. This information is thus at a higher level of abstraction
than the other instance features alone. Contributors to MIPLIB 2017 were invited toprovide complementary material, such as model files in GAMS or AMPL format or
information on block-structured decompositions in the .dec file format. While some
20 submitters accompanied their instances with model files to produce them, no block
structure information was contributed. Nevertheless, GCG, DIP, and DECOMP are
all able to derive information about the block structure automatically (see, e.g., [23,Sec. 5]).
Given the coefficient matrix A of an instance, we formally characterize a (block
structured) decomposition by a partition of the rows Dr = (D1r , . . . , Dkr , L r ) of A
and a partition of the columns D c = (D1c , . . . , Dkc , L c ) of A, such that for every
entry ai, j = 0 of A with i ∈ Dkr 1 and j ∈ Dkc2 it holds that k1 = k2 . We saythat such a decomposition has k blocks. The number k of blocks (normalized by theaverage number of blocks over all instances), the vector statistics (except the standard
deviation) of the numbers of rows and columns per block, respectively, and the socalled area score give rise to ten features related to block-structured decompositions

123

MIPLIB 2017: data-driven compilation...

Fig. 2 Instance b1c1s1 after trivial presolving; left the nonzero entries as they appear in the original
constraint matrix; right with rows and columns re-arranged according to the decomposition produced by
GCG

that were derived from the trivially presolved instances using the detection capabilities
of GCG 3.0.
The area score for a decomposition D = (Dr , D c ) of A is defined askareascore(D) = 1 −

c
r
r
c
r
c
b=1 |Db ||Db | + n|L | + m|L | − |L ||L |mn,which intuitively measures the fraction of the matrix that is “white” in Fig. 2b. A
variant of this score is used by GCG 3.0 to select a decomposition in the case thatseveral, different ones are found (which is typically the case). An area score closer
to 1.0 is better. A very low area score indicates that the model does not contain any
substructures that are easily identifiable by GCG 3.0. This does not imply that there
is no such structure, but rather that it is not obvious.
As computing decompositions can be very time consuming, in particular for huge
instances, GCG was run in a minimal detection configuration. In this configuration,
GCG first groups the constraints of A according to (a) their type in SCIP, (b) their
MIPLIB type (see Table 4), and (c) their number of nonzeros. For each respective
grouping, the smallest cardinality groups of constraints are joined until a specified(small) number of groups remains. From this collection, each subset is potentiallyselected as the set L r of linking constraints of the decomposition. The so-called connected finisher assigns all remaining rows to as many blocks as possible. It is possiblethat this last step identifies a single block only. In this case the area score of the resulting
decomposition is 0.0. From the pool of decompositions that are constructed this way,
one with maximum area score is selected for computing the decomposition featuresof the instance.
4.6 Acquisition of performance data
The selection of the complete MIPLIB 2017 collection (see Sects. 5.2–5.4) was mainly
driven by the feature data described above. The selection of the benchmark set (see123

A. Gleixner, G. Hendel, G. Gamrath et al.
Table 6 List of solvers used for
performance evaluation

SolverVersionThreads

Ref.

CBC

2.9.8

4

[12]

IBM CPLEX

12.7.1

4[13]

Gurobi

7.5.1

4

[25]

MATLAB

R2017b

1

[36]

MOSEK

8.1.0.30

4[41]

SAS/OR

14.2

4

[44]

SCIP

4.0.0

1

[45]

FICO Xpress

8.2

4

[51]

Sects. 5.1, 5.5, and 5.6), however, took into account information about the computational and numerical difficulty of the instances. In order to quantify these empirically,
we collected performance data using current solver software.
Every submitted instance was processed with each of the eight solvers listed inTable 6 to collect performance data. The experiments were performed on two Linuxclusters. The first one consisted of 32 nodes, each equipped with two 3.20 GHz 4-coreIntel Xeon X5672 CPUs and 48 GB RAM, the second consisted of 16 nodes, eachequipped with two 2.50 GHz 10-core Intel Xeon E5-2670 v2 CPUs and 64 GB RAM.
Two such jobs were run in parallel on the same cluster node, each job using a time limitof 4 h and 4 threads, except for SCIP and MATLAB,6 which are both single-threaded.
In total, this performance evaluation required almost 40 CPU years.
Figure 3a shows for every possible cardinality k = 0, 1, . . . , 8 the number of
instances solved by exactly k solvers. 1969 of the 5721 instances (34%) were notsolved by any solver within 4 h. An instance was considered solved if there were
no inconsistencies between solvers and the solution was verified to be feasible (seeSect. 4.7). There were 1155 instances (20%) that could be solved by all eight solvers.
To summarize the results of the experiments, we report here the performance measures for virtual solvers, as described in the introduction. For each instance, a virtual
solver is a summary of all tested solvers by means of an aggregation function such asmin, max, and median, resulting in the best, worst, and median virtual solvers, respectively. The term “virtual” is used to distinguish the presentation from the best (fastest)
or worst (slowest) actual solver over the complete set of instances. The performancemeasures collected are the time to optimality and the number of branch-and-boundnodes processed.
Figure 3b compares the fraction of instances solved by the virtual best and worst
solvers. A large discrepancy between the curves can be observed. The virtual best
solver finished on about 20%, 40%, and 60% of the submissions within 1 s, 1 min,and 1 h, respectively. The virtual worst solver required more than a second for any
instance, and solved only 20.2% of the instances within the time limit of 4 h. The
virtual best solver solved more instances in 2 s than the virtual worst solver was able
to solve in 4 h. Note that all eight tested solvers contributed to the performance of thevirtual best solver, i.e., each solver was the fastest on at least one instance.
6 MATLAB was run using the command intlinprog, which is part of the Optimization Toolbox (TM).

123

MIPLIB 2017: data-driven compilation...
1969

2000

Solved instances

1500

1155
1000

500

385

460

417

348

406
289

291

5

6

0
0

1

2

3

4

7

8

Solvers(a) Number of instances solved by a specific number of solvers.

Frac. of Instances

0.6

virtual solver

0.4

best
worst

0.2

0.0
1 sec

1 min

1h

2h 3h4hTime(b) Fraction of instances solved by virtual best and worst solver for different time limits.
1.00

Frac. of Instances

0.75

Virtual solver (nodes)
best

0.50

worst

0.25

0.00
1e+00

1e+01

1e+02

1e+03

1e+04

1e+05

1e+06

1e+07

1e+08

Nodes

(c) Explored branch-and-bound nodes for solved instances.

Fig. 3 Aggregated results of the performance evaluation on the entire submission

123

A. Gleixner, G. Hendel, G. Gamrath et al.

Figure 3c summarizes data about the number of branch-and-bound nodes processed.
As expected, the number of branch-and-bound nodes varies significantly between
instances, but also between solvers on individual instances. In Fig. 3c, the minimum
and maximum number of explored nodes are shown. In this figure, we consider only
runs that completed sucessfully. Note that there are differences in how solvers report therequired number of branch-and-bound nodes. Concretely, the solution of an instanceduring presolving may be reported as 0 or 1 nodes depending on the solver used. We
therefore normalize the node results so they are always at least 1, i.e., we consider
presolving as part of the root node solution process. This is also justified because weonly consider instances that could not be solved completely by trivial presolving. At
the left end of the scale are the instances that could be solved within 1 node. This
group amounts to 1507 instances, which corresponds to 40% of the instances that
could be solved at all and 26% overall. For more than 50% of the solved instances, thesolution process required less than 1000 nodes. The maximum number of explorednodes is considerably larger. Less than 25% of the considered instances were solved
within 1000 nodes by all solvers that finished within the time limit. Note that for all
385 records (see Fig. 3a) for which only one solver finished within the time limit, theminimum and maximum number of explored nodes coincide.
4.7 Consistency check of solver results
In order to identify numerically challenging instances and incorrect answers returnedby the solvers, the results of the performance runs were independently verified in two
different ways.
First, the feasibility of every primal solution reported at the end of a solution process
was checked. To accomplish this step, a solution checker, whose purpose is to validate
feasibility of the solutions computed by a given MIP solver against the original model,
has been bundled with the MIPLIB scripts since MIPLIB 2010. The checker tries torecognize incorrect results, while at the same time taking into account that most MIP
solvers use floating-point arithmetic, and thus exact feasibility cannot be expected. The
overall structure of the solution checker has been largely unchanged since MIPLIB
2010. It is important to emphasize that it only checks feasibility of a given solution, it
cannot check optimality in general. The solution vector returned by the solver is verified
against the original instance in MPS format, with all computations performed using
the arbitrary precision arithmetic package GMP [24]. Feasibility of linear constraints,
integrality and the objective value are all verified according to given tolerances. We
refer to the MIPLIB 2010 paper [32] for more details on the solution checker designand limitations.
For MIPLIB 2017, we updated the solution checker such that it uses a more flexible(and forgiving) definition of violation of a linear constraint. The previous version of thechecker used an absolute tolerance of ε = 10−4 , so that when given a linear constraintaix ≤ u iA and a solution x ∗ , it would have considered the constraint satisfied if andonly ifaix ∗ − u iA ≤ ε.

123

MIPLIB 2017: data-driven compilation...

However, this turns out to be too strict if the linear constraint has coefficients with alarge absolute value. At the same time, switching to a purely relative tolerance was notconsidered a viable option, as it can lead to a too small effective tolerance when tiny
coefficients are involved. So, in the new version we introduced a different definitionof violation, that tries to combine the strengths of absolute and relative tolerances,
and also to possibly cope with cancellation effects when evaluating the variable part
(activity) of the constraint. In particular, we split the value aix ∗ into its positive
and negative parts as aix ∗ = (aix ∗ )+ − (aix ∗ )− and consider the linear constraint
satisfied if and only ifaix ∗ − u iA ≤ ε · max{(aix ∗ )+ , (aix ∗ )− , |u iA |, 1}.

Given the relaxed definition of violation, it was decided to use the stricter value ε =10−5 as tolerance. The very same logic is applied when checking variable bounds,
and when checking the objective value c∗ reported by the solver—we just treat it as
the linear constraint cx ∗ = c∗ . Note that integrality is still checked with a purely
absolute tolerance of 10−4 . Finally, the solution checker was extended to support
indicator constraints.
Following the feasibility check, which can be done independently for each solver
and solved instance, the results were compared between solvers to identify discrepancies in the optimal values reported or inconsistencies in primal and dual bounds. We
used the publicly available tool IPET [27] to parse the solver log files and validate theresults. We considered results inconsistent when solver A reported a verified, feasible
solution with value c∗ , while solver B timed out reporting a dual (lower) bound that was
higher than c∗ . This included the special case that an instance had been reported infeasible by solver B. For example, on the instance bc1, seven of eight solvers agreed
on an optimal value of 3.338 after exploring search trees with 3k–20k nodes. The
eighth solver, however, reported a solution of value 3.418 as optimal after 720 nodes.
The eighth solver cut off the optimal solution. Note that while such a behavior can becaused by a bug in the solver, it is also possible that different optimal values can be“correctly” obtained when different tolerances are used. Since all MIP solvers rely on
floating-point arithmetics and use feasibility tolerances, the definition of “the optimal
objective value” for a problem instance is ambiguous. In particular, for numericallychallenging problems, a solver might return a different optimal objective value as aresult of applying slightly stricter tolerances within the algorithm. Instances exhibiting such ambiguity are not suitable for benchmarking, since handling this numericalambiguity can be done in different ways, requiring different amounts of computational
effort. This leads to difficulties in comparison. Therefore, we disregarded all instanceswith such inconsistencies during the selection of the benchmark set (see Sect. 5.6),
unless the inconsistency was obviously caused by a bug in one solver; 328 instances(5%) were removed for this reason.

123

A. Gleixner, G. Hendel, G. Gamrath et al.

5 Selection methodology
Due to the vast number of collected instances and the stark overrepresentation of someproblem classes and instance types, it was crucial to reduce the submitted instances
to a carefully chosen selection that provides both researchers and practitioners witha meaningful basis for experimental comparison. MIPLIB 2017 provides two main
instance sets, namely the benchmark set and the collection. In the following, we discuss
the actual selection process and the obtained result.
We approached this task in reverse order by first selecting the larger MIPLIB 2017
collection from the submitted instances, and then choosing the MIPLIB 2017 benchmark set as a subset of the collection. An overarching goal for both the collection andbenchmark sets was to provide good coverage of the feature space of all submissions
while maintaining balance. Note that we thus explicitly avoid allowing the distribution of instance properties observed in the set of submitted instances to affect the
distribution in the final collection, since it is to be expected that the set of submitted
instances would be highly unbalanced in its instance feature representation, if for noother reason than that some submissions contained many more related instances thanothers. A second overarching goal was to choose a collection as large as possible, inorder to obtain a rich and diverse test set, but without sacrificing balance. As explained
above, the benchmark selection step required additional restrictions. With respect tothe benchmark set, the goal was to choose a large set of instances, but with a bias
towards instances that are currently hard for all solvers and keeping in mind that itshould be possible to perform benchmarking in a “reasonable” amount of time.
It seems quite natural to formulate the selection task as an optimization problem.
In fact, we approach the generation of MIPLIB 2017 with a sequence of optimization
problems: a set of diversity preselection MIPs, the collection MIP, and finally, thebenchmark MIP. After an initial cleanup, large model groups (see Sect. 3.4) are cut
down to a handful of diverse instances by the application of the diversity preselection model described in Sect. 5.2. The main purpose of the first selection procedure
was to avoid overrepresentation of instance types from large and very homogeneous
submissions that do not add to the diversity of the instance library. Sections 5.3 and5.5 introduce the clustering procedures to partition the instances based on instancefeatures and performance data, respectively. Sections 5.4 and 5.6 describe the mixed
integer optimization models used to compute the MIPLIB 2017 collection and benchmark sets. Although the selection of the benchmark set is only the final step of theprocess, the initial reduction steps must be aware of which instances are judged tobe suitable for the benchmark set. Otherwise, too many benchmark-suitable instancesmight be excluded initially for selecting a good benchmark set later. Hence, we startin Sect. 5.1 by giving our definition of benchmark suitability.
Note that before the selection steps outlined here, the submission pool of 5721
instances was already reduced to 5666 instances by the removal of LPs (no discretevariables) and instances that are empty after a trivial presolving (see Sect. 4.1). Furthermore, we removed pairs of duplicate instances identified during the selection process.

123

MIPLIB 2017: data-driven compilation...

5.1 Benchmark suitability
The following definition characterizes the requirements for an instance to be in thebenchmark set of MIPLIB 2017. It is important to point out that because we allow
infeasible instances, successful “solution” of an instance for the purposes of this definition means that the solver either produced a (claimed) optimal solution or a (claimed)
proof of infeasibility.
Definition 2 (Benchmark-suitable instance) We call an instance i ∈ I benchmarksuitable if(B.1) it can be solved by at least one considered solver within 4 h;(B.2) it requires at least 10 s with 50% of the solvers;(B.3) it has a constraint and objective dynamism of at most 106 (see Sect. 4.3);(B.4) the absolute value of each matrix coefficient is smaller than 1010 ;(B.5) the results of all solvers on i are consistent (see Sect. 4.7);(B.6) it has no indicator constraints (see Sect. 4.4);(B.7) it is either infeasible or has a finite optimum;(B.8) it is either infeasible, or its (objective) solution value is smaller than 1010 ;(B.9) it has at most 106 nonzero entries.
The subset of benchmark-suitable instances from the ground set I sub is denoted by B.
(B.2) eliminates instances from the benchmark selection that are too easy. Conversely, (B.1) ensures that benchmark instances can be solved by at least one solver, as
already done for MIPLIB 2010. This avoids the situation of MIPLIB 2003, for whichfour instances still remain unsolved 15 years after the release of the test set. The criteria(B.3), (B.4), (B.5), (B.8) ensure that the benchmark set does not contain numerically difficult instances for which results may be ambiguous. Furthermore, benchmark
instances should not contain special constructs that are not supported by all solvers. As
noted in Sect. 4.4, the only special constraint type in the submissions are constraints ofindicator type, which are excluded from the benchmark set via (B.6). (B.7) excludes
feasible instances that do not have a finite optimal value from the benchmark set for
two reasons. First, a feasible, rational MIP has a finite optimal value if and only if
its LP relaxation has a finite optimal value, rendering detection of this property morea continuous than a discrete problem. Second, there is currently no clear consensus
on the expected behavior and output of MIP solvers in the case of a feasible MIP
without a finite optimum. Note that in contrast, infeasible instances are deliberately
not excluded. Finally, (B.9) reduces the hardware requirements to perform tests withthe benchmark set.
Table 7 lists for each criterion the number of excluded instances. Note that an
instance may be excluded for several reasons. In total, 3407 instances were labeled as
not benchmark-suitable, the majority of them because no solver solved them within
the time limit of 4 h.
The larger MIPLIB 2017 collection covers a broader range of MIP instances. It
includes at least one instance from each submitter, a constraint that cannot be enforced
for the benchmark set due to runtime considerations. It may contain instances that123

A. Gleixner, G. Hendel, G. Gamrath et al.
Table 7 Number of instances considered not benchmark-suitable (I sub \ B), as described in Sect. 5.
The column “Crit.” refers to the corresponding criterion in Definition 2. The ground set is the set of 5666
instances available before diversity preselection
Crit.

Exclusion reason

Instances
1958

(B.1)

Too hard: min. solver time > 4 h(B.2)

Too easy: median solver performance ≤ 10 s

741

(B.3)

Objective or constraint dynamism too large

552

(B.4)

Absolute matrix coefficients > 1010

525

(B.6)

Presence of indicator constraints

437
334

(B.5)

Instances excluded for inconsistent results(B.7)

Unbounded instances

87

(B.8)

Best known solution exceeds 1010

80

(B.9)

Too many (> 106 ) nonzeros

40

are considered too easy or too hard for the benchmark set. It may contain instanceswith more dubious numerics suited for testing the robustness of solvers in generaland techniques that explicitly aim at increasing numerical robustness. It may contain
unbounded instances and instances with indicator constraints. It may contain up to
five instances from each model group.
5.2 Diversity preselection
As with previous editions of MIPLIB, the number of instances varies significantly
between different submissions and, more importantly, also between the model groups
described in Sect. 3.4. While some model groups contain a single MIP instance that
represents an optimization problem on a specific data set, other model groups containhundreds of instances using the same underlying model for different data. Hence, forlarger model groups, we preselect a diverse subset of instances as follows.
Let I = I sub denote the index set of submitted instances and let B ⊆ I sub be
the subset of benchmark-suitable instances according to Definition 2. The choice ofa subset of instances can be naturally encoded using a vector of binary variables xi
equal to one if and only if instance i ∈ I is selected. For two instances i, j ∈ I , di, jdenotes the Euclidean distance of their feature vectors. Then for a given model groupG ⊂ I of instances and specified κ ∈ N, we wish to choose κ instances maximallydiverse in the sense that the minimum distance between two selected instances becomesmaximal. If the model group contains benchmark-suitable instances, at least one ofthese should be included in the preselection. Such a preselection can be performed bysolving the mixed binary optimization problem

max

z

s.t.

z ≤ (di, j − d̄)xi x j + d̄

123

(1a)
for all i, j ∈ G , i = j

(1b)

MIPLIB 2017: data-driven compilation...

Diversity Preselection

MIPLIB 2017 Collection

MIPLIB 2017 Benchmark Set

Fig. 4 The results of diversity preselection for the model group “drayage”, which contains a total of 165
instances, in a t-SNE plot. A red diamond indicates that the instance has been selected for the correspondingset

xi = κ

(1c)i∈G

xi ≥ 1

if G ∩ B = ∅(1d)i∈G ∩Bx ∈ {0, 1}G , z ∈ [0, d̄](1e)

An instance i ∈ G is preselected if and only if the corresponding binary variable xi
equals one. The value d̄ := max{di, j : i, j ∈ G } acts as big-M in Constraint (1b). In
order to solve this optimization problem with a MIP solver, the bilinear product xi x j
in (1b) must be linearized by replacing it with an auxiliary binary variable wi, j under
the additional constraints wi, j ≤ xi , wi, j ≤ x j , and wi, j ≥ xi + x j − 1 for all pairsi = j ∈ G .
This preselection was performed for each model group exceeding five instances.
The value of κ used was κG = 5 for groups with 6 ≤ |G | ≤ 10 and κG = 10 for larger
groups. The number of variables and constraints of the preselection model depends
on the size of the corresponding model group. For the largest model group “cmflsp”,
which comprises 360 instances of a capacitated multi-family lot-sizing problem, thediversity preselection instance has 129,242 rows, 64,981 columns, and 323,485 nonzeros. Except for this largest model group, which required approx. 800 s, all preselection
problems could be solved to optimality within a time limit of 500 s using Gurobi 7.5.1.
As an example, Fig. 4 depicts the results of the preselection procedure for themodel group “drayage”, which consists of 165 instances in total. The plot shows
the instances from this group three times. The x and y-coordinates are computed
using t-SNE [50], a technique to embed points of the high-dimensional feature spacein 2D based on their distances. The leftmost plot highlights the optimal solution of123

A. Gleixner, G. Hendel, G. Gamrath et al.

the corresponding diversity preselection instance. Visually, the selected solution forthis group is scattered evenly across the feature space of this group. The middle
and right plot show the instances from this group that are selected for the collection and benchmark set, respectively, for which stricter cardinality restrictions onmodel groups apply. Despite those cardinality restrictions, the corresponding selections appear evenly spread across the ground set. A look at the performance results,
which are taken into account for the selection of the benchmark set, reveals that thetwo selected instances also vary substantially regarding solution time. The easier of
the two instances, namely drayage-100-23 could be solved by seven solvers witha median running time of 40 s. The harder instance drayage-25-23, on the other
hand, could only be solved by three solvers and hence has a median running time of
4 h. The three other instances that are part of the collection lie in between. In total,
diversity preselection reduces the instance set from 5666 to 2182 instances.
The preselected instances form a reduced index set I pre , which serves as input
for the selection of the MIPLIB 2017 collection. Although the following selection
procedure could have been applied to the entire set of submissions, we noticed several
benefits of preselection empirically. It improves the results of the k-means clustering
heuristic described in Sect. 5.3, reduces the size and difficulty of the selection MIPs
to be solved, and finally leads to a larger collection and benchmark set.
5.3 Preparing multiple clusterings
One major challenge in selecting a test set is how to navigate the trade-off of goodcoverage of all observed instance properties against a balanced selection that avoids
overrepresentation. The first goal is best achieved by simply selecting all suitableinstances, while balancedness explicitly asks for removing instances from overrepresented problem classes.
A straightforward method would be to compute one clustering according to theentire feature matrix and pick instances uniformly from each cluster. When applied ina high-dimensional feature space, as in our setting, this naïve approach suffers fromseveral problems well-known in data analysis, such as the curse of dimensionality [5]and the relative scaling of numerical features. The first term refers to the fact that withincreasing dimensionality of the feature space, the difference of distances of one point
to its nearest and to its farthest neighbor becomes smaller in relative terms. Hence,
similar instances cannot be identified reliably. Conversely, depending on scaling, thedistance with respect to one crucial feature may be dominated by less useful features,such that different instances cannot be distinguished reliably. Arguably, the sameproblem holds true everywhere where we use Euclidean distances in the selection
process, for example during the diversity preselection in the previous section. An
important difference between this section and the preselection is that for preselection,
we only considered one model group (with many instances) at a time so that we alreadyknew that all instances were similar to each other. In this and the following sections,however, this model group association is no longer used for the clustering, as we now
incorporate also instances with no known model group association.

123

MIPLIB 2017: data-driven compilation...
Table 8 Clustering statistics for the collection MIP in Sect. 5.4
Feature group

Quality [%]

Lk

Small

Large

δ k,p
Min.

Max.

Variable bounds

91.28

23

7

1

6.57

22.03

Matrix coefficients

87.73

41

2

1

5.33

13.95

Matrix nonzeros

89.10

33

3

0

5.86

16.70

Decomposition

99.99

9

1

1

11.54

15.53

Row dynamism

91.70

17

3

3

9.85

18.12

Constraint classification

87.77

35

6

1

6.63

12.92

Objective nonzero density

93.17

9

0

1

9.95

14.01

Objective coefficients

96.51

43

2

3

4.39

23.05

Sides

98.33

35

10

0

1.00

18.77

Size

87.46

9

2

0

10.95

15.01

Variable types

95.26

7

0

1

8.89

13.54

We therefore counteract the problems of a high-dimensional feature space by using
multiple clusterings of the entire preselected instance set I pre according to disjoint
groups of features. Subsequently, we select instances such that they are balanced withrespect to each of these clusterings. This selection process is more complex and cannot
be achieved by simply picking uniformly from each cluster of each clustering. Instead,
we formulate a mixed integer optimization problem with approximate balancing constraints for each of the multiple clusterings.
Formally, for a given index set I of instances we have K different clusterings, i.e.,
I = Ck,1 ∪ · · · ∪ Ck,L k(2)

for k ∈ K = {1, . . . , K }, with disjoint Ck,1 , . . . , Ck,L k being a partition of theindex set I for every k. The number of clusters L k is allowed to vary, since different
subsets of features may require a different number of clusters to achieve a highquality clustering. We denote the index set of all clusters by C := {(k, ) : k =1, . . . , K ,  = 1, . . . , L k }. Furthermore, the cluster sizes contain outliers, which needspecial treatment in order to avoid limiting the size of the resulting test set too much.
Hence, we partition the set of clusters into small, medium (regular-sized), and large
clusters and denote the respective index sets by S, M, and L ⊆ C, as follows. A cluster
Ck,l is denoted small if its size is less than half the average size of Ck,1 , . . . , Ck,L k . On
the other hand, a cluster is treated as large if it is displayed as an outlier in a typical
boxplot. Concretely, Ck,l is considered large if its size exceeds the 75% quantile among
Ck,1 , . . . , Ck,L k by more than 1.5 interquartiles.
For the selection of the MIPLIB 2017 collection, we use one clustering for eachof the K = 11 groups of instance features listed in Table 4. The clusterings of allpreselected instances (2182 instances) are computed using a k-means heuristic [26],which yields a first family of clusterings denoted by K1 = {1, . . . , 11}.

123

A. Gleixner, G. Hendel, G. Gamrath et al.

Table 8 gives insight into the result of this clustering process. It shows the total
number of clusters (value of L k ) for each feature group clustering. The quality column
describes the percentage of the total feature group distance between clusters. Moreformally, for a clustering k ∈ {1, . . . , K } of instances, the quality of this clustering isi= j∈I di, j −L k 
=1i= j∈Ck, di, ji= j∈I di, j· 100 %.

The value range of the above fraction is the interval [0, 1] such that the quality lies
between 0 and 100%. A clustering has a high quality if long distances between
instances are between different clusters. Note that for the distance computation for
the quality measure, only features contained in the corresponding feature group wereconsidered. The table shows that the quality of the clustering was at least 87% andoften significantly above 90%, yielding an average quality of 92%. The individual
value of L k has been manually selected for every feature group as the minimum
integer L that admits a clustering with at least 90% quality over the set I sub (5666instances) unless the targeted quality was not achievable using a reasonable amountof clusters. In addition to the number of clusters L k , the table presents the number ofsmall and large clusters, which are constrained less strictly than the medium clusters,see Constraints (4a) and (4b) below.
The last two columns report the minimum and maximum total dissimilarity perfeature group. For each cluster (k, ) ∈ C, its total dissimilarity δk, ≥ 1 is defined as
the shifted geometric mean of the pairwise Euclidean distances {di, j : i < j ∈ Ck, }between its instances, using a shift of 1. Here, distances are computed with respect tothe entire feature space, in contrast to the above cluster quality computation. Becauseof the shift by 1, the smallest possible value of δk, is 1, which only occurs for clusters
containing exactly one element, or for clusters that contain only instances which are
indistinguishable in the feature space. All cluster parameters from Table 8 enter theselection constraints described in the next section.
For the selection of the benchmark set, we additionally use performance data to
hand-craft three clusterings for each of the eight participating solvers, which yields
additional clusterings K2 = {12, . . . , 35}, see Sect. 5.5 below.
5.4 Selection of the MIPLIB 2017 Collection
In the following, we describe linear formulations to enforce the requirements specified
by the committee. At this stage, the instance set was limited to the instances I = I preleft after the diversity preselection procedure. The set of clusterings K = K1 was
the one determined using the instance feature groups from Table 4.
To express balancedness, consider one clustering I = Ck,1 ∪ · · · ∪ Ck,L k .
Naïvely,
we would like to pick the same number of instances from each cluster, i.e.,
i∈Ck, x i ≈ yk for an auxiliary variable yk ≥ 0. However, enforcing this for all clusterings is highly restrictive. Furthermore, while the instances in each of the clusters
Ck,1 , . . . , Ck,L k should be homogeneous with respect to the features that were used

123

MIPLIB 2017: data-driven compilation...

to compute clustering k, they may be heterogeneous with respect to the entire feature
vector. This interaction between different clusterings must be taken into account.
To achieve this, we consider the total dissimilarity of the clusters that was introduced
above. Arguably,from clusters with higher total dissimilarity, more instances shouldbe picked, i.e., i∈Ck, xi ≈ δk, yk . Introducing a tolerance parameter ε, 0 < ε < 1,
we arrive at the balancedness constraints(1 − ε)δk, yk ≤xi ≤ (1 + ε)δk, yk .

(3)i∈Ck,

Concretely, we used ε = 0.5 for the selection of the collection and the benchmark set.
In practice, we discard the left inequality for small clusters and the right inequality forlarge clusters and use
xi ≥ (1 − ε)δk, yk

for all (k, ) ∈ C \ S,(4a)xi ≤ (1 + ε)δk, yk

for all (k, ) ∈ C \ L.

(4b)i∈Ck,

i∈Ck,

Additionally, if two instances have identical feature vectors, then at most one of themshould be chosen , i.e.,
xi + x j ≤ 1

for all i, j ∈ I × I with i < j, di, j = 0.

(4c)

At most five instances should be selected from each model group. If the modelgroup contains benchmark-suitable instances, at least one of those should be included
into the MIPLIB 2017 collection. Let I = G1 ∪ . . . ∪ G P denote the partition of
instances into different model groups, then this condition readsxi ≤ 5for all p = 1, . . . , P,

(4d)xi ≥ 1for all p = 1, . . . , P with G p ∩ B = ∅.

(4e)

i∈G p

i∈G p ∩B

Furthermore, from each submitter at least one instance should be selected, i.e.,xi ≥ 1for all s = 1, . . . , S,

(4f)

i∈Ss

where I = S1 ∪ · · · ∪ S S denotes the partition of instances with respect to S submitters.
Finally, we imposed relative limits on a small number of specific subsets of123

A. Gleixner, G. Hendel, G. Gamrath et al.
Table 9 Instance sets for which relative limits on the selection apply, with limits shown for the collection
and benchmark MIPs individually. The Size columns show the size of this set within the respective groundset the selection is based on. The parameter ρk is a relative limit on the allowed instances from this set in asolution as specified by Constraint (4g)k

SetsCollection MIP (2182)
Size

ρk

Benchmark MIP (499)
Size

ρk
0.05

1

MiniZinc instances

484

0.05

14

2

NEOS instances

696

0.33

182

1.00

3

small (n ≤ 2000)

691

0.20

124

0.20

4

medium (n ≤ 10,000)

1309

0.50

290

0.50

5

BP (n = n b )

422

0.20

109

0.20

instances, Rr ⊂ I , by requiringxi ≤ ρri∈Rrxi

for all r ∈ {MiniZinc, NEOS, small, medium, BP}.

(4g)

i∈I

The concrete values for ρr , for both the collection MIP and the benchmark MIP
described in Sect. 5.6, are given in Table 9. The numbers in parentheses show the sizeof respective ground sets I pre and I col ∩B, from which instances were selected. The
number of instances from the NEOS server was limited by the committee because of thelack of reliable information on their application and model background. Purely binary
problems (BP) were limited because they often represent academic applications suchas combinatorial puzzles, but less often occur in industrial, “real-world” instances. The
limit ensures that enough actual mixed integer instances are selected for the collection
and benchmark sets. For the groups in Table 9 that refer to instance features, thefeatures are always evaluated after trivial presolving.
Subject to those constraints, our objective was to include as many instances as
possible, preferring benchmark-suitable instances. Hence, we formulate the collection
MIP as the mixed binary optimization problem
βi xi : (4a) − (4g), x ∈ {0, 1}I

max

pre

, y ∈ RK
≥0 ,

(5)

i

with objective cofficients β to prefer benchmark-suitable instances,
βi =2, if i ∈ B,
1, otherwise.

We solved the collection MIP over the ground set I pre of 2182 instances (after diversity preselection). Despite our efforts to remove obvious duplicates, there remained
48 pairs of instances in I pre with a feature distance of zero. From each such pair, at
most one instance was selected for I col because of Constraint (4c). We obtained the123

MIPLIB 2017: data-driven compilation...
Table 10 Influence of the parameter ε used in Constraints (4a) and (4b) on the number of instances selected
for the collection set. The second row reports how many of the selected instances are benchmark-suitable
ε

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Selected instances

–

–

740

946

1065

1110

1125

1140

1140

Benchmark-suitable

–

–

355

465

499

546

560

561

564

MIPLIB 2017 collection I col , which comprises 1,065 instances, 499 of which arebenchmark-suitable.
The choice of the balancedness parameter ε used in Constraints (4a) and (4b) clearlyplays a central role in the formulation of the collection MIP. In order to analyze thesensitivity of the result to this parameter, we solved the collection MIP for different
values of ε. Table 10 reports on the sizes of the resulting collection sets. The smallest
tested value of 0.1 makes the collection MIP infeasible, while the second smallest
value of 0.2 did not return with a feasible selection after 15 minutes. For larger valuesof ε, the collection MIP is feasible and can be solved reasonably quickly. The numberof selected and benchmark suitable instances increases along with the balancedness
parameter and allows us to control the number of selected and benchmark-suitableinstances.
5.5 Performance clusterings
In addition to instance features that depend exclusively on instance data, computationaldifficulty is an important aspect to consider for the benchmark set. We assessed thecomputational difficulty of every instance empirically by considering the performancedata of the eight tested solvers (see Sect. 4.6). To quantify performance, we considereda matrix of running times tw,i > 0 for each solver w ∈ W := {1, . . . , W } andinstance i ∈ I . If w could not solve the instance i, tw,i was set to the time limit of 4h. We denote by Iw ⊆ I the set of instances that were solved by w within the time
limit.
For each of the participating solvers, we created three different clusterings of theinstances to capture different aspects of performance. The base set I for these clusterings are the 499 benchmark-suitable instances within the MIPLIB 2017 collection,i.e., I = I col ∩ B. The overall goal was to avoid a biased selection of instances, i.e.,
to avoid that the absolute and relative performance of a solver on the benchmark set
appears different than on I col . Each of the three clusterings avoids a different bias.

Absolute performance clustering The first clustering uses an absolute performanceranking. For each solver w, we sorted the instances in Iw according to increasing
running time tw,i . For a fixed number of clusters B, which we set to B = 11, wegrouped the instances in Iw into B equally-sized clusters w.r.t. increasing rank, i.e.,
we assigned the instances solved fastest to the first cluster, the next fastest set of
instances to the second cluster, …, and the slowest instances to the last cluster, so thateach cluster contained approximately |Iw |/B instances.

123

A. Gleixner, G. Hendel, G. Gamrath et al.

The instances in I \ Iw that could not be solved by solver w seem indistinguishable with respect to performance. However, it could be that the solution process was
terminated only seconds prior to concluding the proof of optimality, or that the solution process would have continued unfinished for days or even months. We took this
into account by inspecting the performance of the other solvers and formed two moreclusters from those instances: instances that could be solved by exactly one solver andinstances that could be solved by more than one solver. The case that instances could
be solved by no other solver does not appear since such instances are not benchmarksuitable (Definition 2, Criterion (B.1)).
Relative performance clustering In contrast to the absolute performance clustering,the instances were ranked based on relative solver performance for a second clustering
as follows. To this end, we defined the relative performance of solver w on instance iwith respect to the other solvers as
rel
:=
tw,i

tw,i + σ,
min tw ,i + σ

(6)

w =w

where σ ∈ R≥0 is a nonnegative shift as in the computation of shifted geometricmeans. Relative performance locates the individual solver performance relative to allother solvers on an instance, regardless of the absolute scale. The motivation is thatsolvers and solver improvements are traditionally measured by the shifted geometricmean instead of the arithmetic mean. For the instances that could be solved by thissolver, we used this ranking to define B equally-sized clusters in the same fashion aswith the absolute performance ranking. The timeout instances were again divided into
two further clusters of instances that could be solved by exactly one and by more than
one other solver, respectively.
Binned absolute performance clustering The third clustering uses absolute solving
time directly, partitioning possible solving times into B = 7 intervals[T0 = 0, T1 ), [T1 , T2 ), . . . , [TB −1 , TB ),(7)

whose breakpoints are equal for all solvers. The concrete bin width used grows exponentially as follows.
T j = 10−3.5+0.5 j · 14400,

j = 1, . . . , 7.

Hence, the righmost bin T7 has the time limit of 4 h as right breakpoint. Then foreach solver w ∈ W we formed B clusters {i ∈ Iw : tw,i ∈ [T j−1 , T j )}, j =1, . . . , B . Empty clusters were discarded. This is different from the absolute andrelative performance clusterings in that it partitions the instances solved by a solver
into clusters that differ in size. The instances in I \ Iw , which could not be solved
by solver w, were again treated as two further clusters as above.
All in all, this led to 24 clusterings, K2 = {12, . . . , 35} The ranking-based clusterings yield approximately equal cluster sizes over Iw , but these cluster sizes may differ

123

MIPLIB 2017: data-driven compilation...

to the ones on I \ Iw . The binned absolute performance clustering does not control
cluster size and may yield very unequally sized clusters. In the following benchmark
MIP we picked from the performance clusters according to their size, i.e., we use
δk, = |Ck, | for all k ∈ K2 in Constraint (4a) and (4b).
In the case of SCIP, as an example of an absolute performance clustering, each of the
11 parts contained between 20 and 22 instances, and the remaining 263 instances weresplit into 58 instances which could only be solved by one solver, and 205 instancessolved by at least two solvers. Although the absolute and relative performance clusters
were, by design, almost equal in size, we observed quite different partitions of theground set. An example is the fastest relative performance cluster for SCIP, whichshares 8 of its 21 instances with the fastest absolute cluster. The remaining 13 instances,
for which SCIP was particularly fast compared to its competitors, were spread across8 of 10 possible absolute clusters. This shows that, to some extent, the two suggested
clusterings exhibit an almost orthogonal instance partition.
5.6 Selecting the MIPLIB 2017 benchmark set
The benchmark set was selected from the ground set of benchmark-suitable instancesin the MIPLIB 2017 collection, I = I col ∩ B. The balancedness constraints (4a)and (4b) are now defined using instance feature and performance clusterings, K =K1 ∪ K2 . The rationale is that the current performance of solvers in the collection
should be reflected by the performance on the benchmark set in order to avoid any
unintentional bias towards a solver during the selection of instances. The relative limit
constraints (4g) are kept, but the restriction on instances from the same model group is
reduced to one instance from NEOS groups and two instances, otherwise. The stricter
limit on NEOS instances is imposed because little information is available for theseanonymously submitted instances and the chance for duplicate instances in the samemodel group is deemed higher.
In addition, executing one benchmark run on this test set should be possible within areasonable time frame on modern hardware, possibly a compute cluster. We specified
a total time limit τ of 32 days for running the benchmark set with a hypothetical solver
with median running times capped to a time limit of 2 h, i.e., t¯i := min{median{tw,i :w ∈ W }, 7200}. The resulting benchmark MIP reads
βi xi

max

(8a)

i

s. t. (4a), (4b), (4g),

1
xi ≤
2
i∈G p

for all p = 1, . . . , P from NEOS,
otherwise,(8b)t¯i xi ≤ τ,(8c)x ∈ {0, 1}I ,(8d)y ∈ RK
≥0 .

(8e)

i∈I

123

A. Gleixner, G. Hendel, G. Gamrath et al.

This approach has a potential drawback. Representing current solver performance
may overly favor instances that are tractable by current solver technology. This is
opposed to one main goal of the MIPLIB project, which is to provide a test bed thatdrives solver development forward. Hence, we used the objective coefficient
βi := 1 +1min tw,i
20 w∈W(9)for instance i ∈ I . This favors instances that are challenging even for the virtual best
solver. The concrete choice of the square root was empirically motivated.
Using I col ∩ B containing 499 instances as ground set, we solved the benchmark MIP that respects all feature group and performance clusterings. The solution to
the benchmark MIP contained 240 instances and now constitutes the MIPLIB 2017
benchmark set I bench . We note that the imposed running time constraint (8c) was nottight on our performance data. A solver with median running times and a time limitof 2 h would take about 14 days to process all benchmark instances sequentially.
We also note that for several reasons, the goal of representing computational difficulty in the reduced benchmark set cannot be achieved perfectly and is not evenwell-defined: The performance data is only a snapshot of current algorithms. It was
gathered using a time limit, performance variability and parallel scalability were notcaptured, and correctness was only enforced approximately with respect to the tolerance parameter ε. Last but not least, the different performance clusterings may evencontradict each other. However, we hope that it helps to avoid unintentionally andunconsciously introducing a bias both towards instances of a particular difficulty andtowards any of the solvers.

6 The ﬁnal collection and benchmark sets
The MIPLIB 2017 collection I col has been compiled with a focus on a balanced
and diverse representation/coverage of the feature space. The benchmark set I bench
incorporates similar requirements also for the performance data. This section discusses
the feature and performance aspects of the compiled sets. We also assess the descriptive
power of the feature space by (re-)detecting known model group associations.
6.1 Representation in feature space
A frequent question during the discussions about the MIPLIB 2017 compilation process was whether the choice of features and their scaling is able to distinguish instances
in a useful way. Ideally, two instances based on the same model for the same application, but with different data, should be close to each other in the feature space,regardless of variations of, e.g., the size of the matrix. For MIPLIB 2017, we have twosources to assess similarity between instances, namely their distances in feature space
and the model groups G from Sect. 3.4. In this paragraph, we evaluate the descriptive
power of the feature space by comparing similarity in feature space and model groupassociation of instances.

123

MIPLIB 2017: data-driven compilation...

Let X denote the instances in the MIPLIB 2017 collection (|X | = 1065). The
complete graph K X = (X , E) on the vertex set X has |E| = |X2 | = 566, 580 edges.
Now consider two subsets of the edges. Let E G denote edges between instances fromthe same model group. Furthermore, let for every x ∈ X , Sx ⊂ X \ {x} denote theset of five most similar instances to x in the collection w.r.t. the distance in the scaledfeature space. With the sets Sx , we define E S to be the set of similarity edges as
E S := {(x, y) ∈ E : x ∈ S y or y ∈ Sx }.
Note that x ∈ S y does not necessarily imply the opposite containment y ∈ Sx , butholds in many cases. The actual cardinalities of the two sets are |E G | = 1327 and|E S | = 3747 and hence comprise less than 1 % of the total edge set. Indeed, computing
the probability for an edge e that was selected uniformly at random from E to be asimilarity edge isP(e ∈ E S ) =|E S |≈ 0.007|E|

Now what is the probability that a group edge is also a similarity edge? This question
can be answered by computing the conditional probabilityP(e ∈ E S | e ∈ E G ) =|E S ∩ E G |974
P(e ∈ E S ∩ E G )==≈ 0.734.
P(e ∈ E G )|E G |
1327

The majority of group edges is contained in the similarity set, and the probability for agroup edge to be a similarity edge is more than 100 times higher than for a randomly
selected edge.
Recall from Sect. 3.4 that the model groups have been partially derived from thefeature data. For submissions from the NEOS server, which have an unknown origin,clustering has been used to group similar NEOS instances into pseudogroups. If weomit all NEOS instances from the above computations (by considering the completegraph K X \X NEOS with 714 vertices), the probability for an edge to be a similarity edge
is about the same, P(e ∈ E S ) ≈ 0.008, whereas the conditional probability of a groupedge to be a similarity edge is ≈ 0.815 and hence even higher than for K X .
From this observation, we conclude that the feature space has been designed sufficiently well for the clustering approach. In fact, the used feature space recovers themodel group data better than we expected. Even for an instance that does not belong toa dedicated model group or lacks bibliographical information, the similarity to othermodel groups can yield interesting hints at the type and application of this instance.
Therefore, the web page of MIPLIB 2017 (see also Sect. 6.3) allows to browse thefive most similar instances for every instance of the MIPLIB 2017 collection.
Figure 5 uses t-SNE [50] to give a spatial impression of the locations of the MIPLIB
2017 benchmark set, the benchmark-suitable instances, and the collection, relative toeach other in feature space. The distance computation is based on the feature vectors
after they have been scaled over the entire set of submissions. Note that there is a123

A. Gleixner, G. Hendel, G. Gamrath et al.

MIPLIB 2017
Benchmark Set (240)Benchmark−suitable (499)Collection (1065)

Fig. 5 Two-dimensional embedding of the distances in the feature space using t-SNE

subset relation for those sets, i.e.,
Benchmark Set  Benchmark-suitable  Collection,such that the plot only shows the innermost set membership for every instance. Figure 5
shows that benchmark-suitable instances cover the majority of the feature space thatis comprised by the collection except for a few regions.
6.2 Solver performance
One of the goals of MIPLIB has always been to provide a representative set to measure
and compare solver performance. In this section, we analyze the solver performanceon the collection, and in particular on the new benchmark set of 240 instances. For
this analysis, we use the computational results obtained during the 4 h runs conducted
for the selection process. Note, however, that in this article, solver performance is notreported directly for several reasons. One reason is that due to hardware restrictions,
not all runs could be performed exclusively on the same hardware, and should hence
not be reported in a way that could be confused for an actual benchmark. Again, theindividual results are aggregated into the virtual best solver, i.e., a solver that always
finishes as fast as the fastest actual solver for each individual problem instance.
In Fig. 6, we compare the performance of this virtual best solver on the benchmark
sets of MIPLIB 2010 and 2017. One of the motivations for the creation of MIPLIB
2017 was the demand for a harder benchmark set. As a consequence of Definition 2,
the virtual best solver solves all instances within 4 h (or 14,400 s) as this is one of
the criteria for benchmark suitability. The plot shows that the majority of the oldbenchmark set can be solved by the virtual best solver in less than 100 s, and that there
is no instance left where the virtual best solver requires 1 h or more. By contrast, thebenchmark set of MIPLIB 2017 is much more demanding, as a significant portion of123

MIPLIB 2017: data-driven compilation...

15

10

Instances

Set
MIPLIB 2010 (87)
MIPLIB 2017 (240)
5

0
1

10

100

1000

10000

Min Time (sec.)

Fig. 6 Minimum solving time in seconds of the benchmark sets of MIPLIB 2010 and MIPLIB 2017

Table 11 Percentage ofbenchmark instances that could
not be solved within 1–4 h by thevirtual best and median solvers

Virtual median2010 (%)

Virtual best2017 (%)

2010 (%)

2017 (%)
19.6

>1h

17.2

66.7

0

>2h

12.6

61.7

0

8.8

>3h

6.9

51.2

0

4.2

>4h

5.7

48.3

0

0.0

The percentages are relative to the individual benchmark sets (2010:
87 instances, 2017: 240 instances)instances lies at the right end of the scale. Due to its increased size, the bars for the
MIPLIB 2017 benchmark set lie almost consistently above the ones for the previous
set. The MIPLIB 2017 benchmark set covers much more of the relevant performancescale than its predecessor covers nowadays. Note that the MIPLIB 2010 benchmark set
appearing easy is an impressive result of 7 years of continuous solver improvements.
Table 11 shows the percentage of instances of the respective benchmark set (2010or 2017) for which the virtual best solver takes longer than a certain time threshold,which we vary between 1 and 4 h. As mentioned before, all instances of the 2010
benchmark set could be solved within 1 h by the virtual best solver. The table showsthat among the new benchmark set, almost 20% of the instances cannot be solved
within 1 h by the virtual best solver. Along with the virtual best solver, Table 11 also
shows the performance of the virtual median solver, based on the median solution
time over the eight involved solvers. Since the number of involved solvers is even, themedian is computed by averaging the timing result of the two solvers ranking 4th and5th for each instance individually. Therefore, the virtual median solver is faster than123

A. Gleixner, G. Hendel, G. Gamrath et al.

half of the solvers. On the MIPLIB 2010 benchmark set, 17.2% of the instances are
not solved within 1 h by the virtual median solver. In contrast to the virtual best solver,
the virtual median solver still times out after 4 h on 5.7% (5 of 87) instances evenon the 2010 benchmark set. On the MIPLIB 2017 benchmark set, the virtual mediansolver needs more than 1 h on two thirds of the instances, and still needs more than 4
h of solving time on almost 50% of the instances.
In Fig. 7, the fraction of instances solved by the virtual median solver as a functionof time is shown. The figure shows the corresponding curves for the MIPLIB 2017
collection of 1065 instances, the set I col ∩ B of 499 benchmark-suitable instances,and the MIPLIB 2017 benchmark set (240 instances). Recall that an instance is only acandidate for the benchmark set if it takes the virtual median solver at least 10 s to solve
it, which is also visible from the plot. As minimum for any time measurement, 0.5 s
were used, which is visible in the curve of the MIPLIB 2017 collection. The objectivefunction for the benchmark set was designed to prefer harder instances. The effect of
this design choice is visible in the figure, in which the percentage of solved instancesof the benchmark set consistently lies below the curve for the 499 benchmark-suitableinstances. The slope of the curve for the collection is approximately linear for about
90% of the visible area. Due to the logarithmic scaling of the horizontal (time) axis,
this suggests that to a certain extent, the solving behavior of the virtual median solver
can be approximated by fitting a logarithmic function. Note that the clear change inbehavior of the curves, which are otherwise almost logarithmic, towards the right endof the scale is an artifact from the median computation, which weighs in as soon as the
4th solver could still solve the instance, but the 5th solver couldn’t. Their timings caneven be very different. On the instance blp-ic98, the four best performing solvers
finish within 1700 s, while the fifth solver times out after 4 h. The median solver
therefore has a performance of 8050 s. All individual curves of the actual solvers
tested have a similar, almost logarithmic shape without the median artifact.
Statistically speaking, the curves in Fig. 7 describe empirical cumulative density
functions (CDF) of the random variable that represents median solving time for an
instance. The Kolmogorov–Smirnov (KS) test is a statistical approach to measuring
the similarity between a pair of CDF F1 , F2 . To this end, the KS test measures the
maximum vertical distance.7 D between F1 and F2 . With increasing D, the likelihood
decreases that F1 and F2 represent samples from the same distribution. In order tofurther quantify the hardness of the benchmark set, a KS test has been performed forevery solver, comparing its individual CDF pair on the benchmark set and the set ofbenchmark-suitable instances. As alternative serves the hypothesis that the CDF of
the benchmark set lies below the CDF of the larger set.
For the performance of the virtual median solver as depicted in Fig. 7, the distance
D is approximately 0.10, which results in a p value of 0.04. A much smaller p-valueof 6.513 · 10−5 is obtained for the virtual best solver at a distance of D = 0.17. For
four of the actual solvers, the maximum D is larger than 0.1, and the accompanying
p value is smaller than 1%. This is significant evidence that for those solvers, thebenchmark set has been selected as a particularly hard selection among all suitableinstances. For the other four solvers, the value of D is smaller, which in turn results
7 This distance D is the supremum norm sup

123x∈R |F1 (x) − F2 (x)|.

MIPLIB 2017: data-driven compilation...

1.00

Frac. of Instances

0.75

Set
Benchmark−suitable (499)0.50

MIPLIB 2017 Benchmark Set (240)
MIPLIB 2017 Collection (1065)

0.25

0.00
10

1000

median Time (sec.)

Fig. 7 Virtual median solver performance on subsets of the MIPLIB 2017 collection

in larger p-values (each greater than 0.1). Note that even here, the curves of the CDF
on the benchmark set tends to undercut the CDF of the set of suitable instances, but
this discrepancy is not large enough to render the KS test significant. Hence, for thosefour solvers, the performance curve is more or less representative of the CDF overall suitable instances. The results show that the selection methodology has achieved
both its conflicting goals, hardness and representability, with respect to the solver
performance, equally well.
6.3 The MIPLIB 2017 Web Page
For the release of MIPLIB 2017, the web page http://miplib.zib.de has been writtenfrom scratch and received a modern design. The main page shows the current status of
instances regarding the categories easy, hard, and open. The two main tables list the
instances of the MIPLIB 2017 collection and benchmark set together with some key
properties, their model group, and their optimal or best known objective value. All
tables use tags on the instances to highlight certain properties that may be interesting for
researchers, such as pure feasibility instances with no objective function, instances forwhich good decompositions are known, instances with critical numerics, the presenceof indicator constraints, etc. It is possible to search and filter for instances by name,status, model group, or tag, or to sort the table by column.
The individual instance pages offer a short description of each instance and bibliographical information. Also, more information on the constraint mix of each instance
before and after presolving is displayed, as well as decomposition information, if available. Finally, the optimal or best known solutions for every instance are displayed, as

123

A. Gleixner, G. Hendel, G. Gamrath et al.

well as the five most similar instances as explained in Sect. 6.1. The web page alsooffers additional downloadable content, including– the MIPLIB 2017 collection and benchmark sets,– lists of instances with certain tags,– available solutions,– optimal/best known objective values for all instances,– the collection and benchmark MIPs,– the feature extractor, and– bash scripts to run and validate solver performance experiments.
Some of the tags may change over time, such as the instances that fall into thecategories easy, medium, and hard. Also, best known objective values are naturally
changing or eventually proven to be optimal. Therefore, we provide versioned files for
accurate referencing. The versioned files are periodically updated. All downloadable
solutions are checked for feasibility with the solution checker from Sect. 4.7. Also,
their exact solution values after fixing all integer variables are computed using SoPlex
with iterative refinement. The actual collection and benchmark MIPs are also available
for download. Obviously, these instances cannot be part of the actual collection andbenchmark sets, respectively, since their presence would alter the feature space andhence their own formulation. Contributions in terms of updated bibliographic information or corrections to the instance descriptions are very welcome. In particular, weare constantly accepting and checking improving solutions to the open instances ofthe MIPLIB 2017 collection. In contrast to previous MIPLIBs, not only new optimal,
but any improving solution will be considered for the periodic update of the page data.
Improving solutions should be sent to the maintainers of the page, together with adescription of how they have been obtained. Note that every submitted solution must
adhere to the format accepted by the MIPLIB solution checker (Sect. 4.7), which isalso available on the web page.

7 Conclusion
The sixth version of MIPLIB has, as was the case in previous updates, significantly
increased in size compared to its predecessors. The distinction between a dedicatedbenchmark set and the entire collection, which was introduced with MIPLIB 2010,
has been preserved. These sets now contain 240 and 1065 instances, respectively.
The process of how to reduce the initial submission pool of over 5000 instances to abalanced selection of this size, however, has been completely redesigned. Beyond the
new MIPLIB 2017 itself, the development of this fully data-driven and optimizationbased methodology is the main contribution of this paper.
We propose two related MIP models that have successfully provided decision support for the selection process to the MIPLIB committee. One key ingredient of thisapproach is the definition of a feature space covering more than a hundred different
dimensions for characterizing a MIP instance. In order to ensure a balanced selection
with respect to these features and, for the benchmark set, with respect to performancedata, we advocate the use of multiple instance clusterings over partitions of the fea-

123

MIPLIB 2017: data-driven compilation...

ture vector. A comparison with manually assigned model groups available from metadata of the submissions shows the high descriptive power of the used feature space.
By approaching the selection problem as a MIP that encodes a balanced, simultaneous selection from each such clustering, the collection and benchmark MIPs provide
the flexibility to incorporate both feature coverage and performance requirements aswell as other restrictions from the comittee. Besides improving the final outcome, thisformalization of the selection criteria has served to increase the transparency of theselection process.
While the selection methodology proposed here is not intended as a general
blueprint for construction of test sets, we hope that parts of the process of constructing
the MIPLIB instance sets may apply to the curation of other test sets in the future.
Certainly, many variations and adjustments of the approach are possible. Not only
the chosen constants or heuristic clustering methods can be adapted, but also the roleof objective function and constraints may be redefined. Furthermore, the interplay
between the main selection models and the diversity preselection offers potential forvariation. For example, a different approach may directly select and fix a number of
instances with maximum diversity from each large model group for the collection andafterwards complete the collection with instances from smaller groups and instanceswith no known model group association. In light of these degrees of freedom and many
ad hoc decisions that had to be made in advance, the final result is clearly only oneof many possible and justified outcomes. However, we believe that the collection andbenchmark sets presented in this paper are a profound attempt to provide the researchcommunity with test sets that represent the versatility of MIP.
One of the main characteristics of the benchmark set is that it provides a snapshotof current solver performance. Our hope is that the performance of solvers on thisbenchmark set will serve as a sort of bellwether for progress in the field of mixed
integer optimization as a whole in the coming years. As future work, we propose
to assess the performance representability of the benchmark set from hindsight, i.e.,
by comparing speed-ups for entire model groups as well as for the selected instances.
Such data will finally allow to better compare different (pre-)selection models such as,e.g., the one presented here that favors diversity to a different one that selects nearest
neighbors and maximizes representability.
Another benefit of our MIP-based selection process is the fact that the MIP models can be used to approach questions beyond the initial creation of the test set.
One example is the following case, in which it occurred that benchmark instances
needed to be replaced as new computational data became available. Despite all thecare that was taken to exclude numerically critical instances from the benchmark set,problematic numerical properties of the two instances neos-5075914-elvire
and neos-3754224-navua remained undetected during the selection process.8 A
variant of the benchmark selection MIP was used to compute a minimal update of
the benchmark set that exchanges the discussed instances while preserving the bal8 Both instances are at the border between feasibility and infeasibility, but at the time of collecting solution
data no inconsistencies could be observed. For the first instance, two solvers agree on the optimal solution
value although the instance should mathematically be infeasible. The second instance has only been declaredinfeasible by one solver during the selectiong process; we received a solution that is feasible within tolerances
half a year after the original publication of the benchmark set.

123

A. Gleixner, G. Hendel, G. Gamrath et al.

ancedness requirements. An accordingly updated version of the benchmark set was
published in June 2019.
Finally, by the time of this writing, the challenges of MIPLIB 2017 collection havealready attracted a wide audience. In fact, we have received many new solutions topreviously open instances. While some of those optimality or infeasibility proofs havebeen obtained by the use of massively parallel codes such as the Ubiquity Generator
framework [46,47], other instances inspired the development of customized cutting
plane approaches, or even a rigid mathematical proof of infeasibility without any codein the case of the instance fhnw-sq3. In total, 68 originally open instances havealready been solved.9 We are looking forward to further contributions and many more
years (and versions) of MIPLIB to come.
Acknowledgements The authors wish to thank all members of the MIPLIB 2017 committee, the creatorsof the previous editions of MIPLIB, all submitters to previous MIPLIB collections, and most importantly
all contributors of new instances for their effort: Alexandra M. Newman, Andrea Arias, Andreas Bärmann,
Andrew Stamps, Antonio Frangioni, Austin Buchanan, Balabhaskar Balasundaram, Berk Ustun, Cezar
Augusto Nascimento e Silva, Christian Liebchen, Christopher Daniel Richards, Christopher Hojny, Dan
Hiroshige, Koichi Fujii, Dan Neiman, Daniel Bienstock, Daniel Heinlein, Daniel Rehfeldt, Dimitri Papageorgiou, Domenico Salvagnin, Felix Cebulla, Felix J.L. Willamowski, Gavin Goodall, George Fonseca,
Gerald Gamrath, Gerald Lach, Gleb Belov, Hans Mittelmann, Haroldo Gambini Santos, Hsiang-Yun Wu, Irv
Lustig, Janos Hoener, Jeff Linderoth, Jesus Rodriguez, Jonathan Eckstein, Jordi Castro, Joshua Friedman,
Juan Javier Dominguez Moreno, Koichi Fujii, Laurent Sorber, Manuel Iori, Marc Pfetsch, Mark Husted,
Matias Soerensen, Michael Bastubbe, Michael Winkler, Paula Carroll, Pelin Damci-Kurt, Philipp Leise,
Pierre Le Bodic, Qie He, Rob Pratt, Salim Haddadi, Sascha Kurz, Sean MacDermant, Shunji Umetani, Simon
Bowly, Simon Felix, Siwei Sun, Stephan Beyer, Sujayandra Vaddagiri, Tamas Terlaky, Timo Berthold, Toni
Sorrell, Utz-Uwe Haus, and Yoshihiro Kanno. Furthermore, we are grateful to the three reviewers and thearea editor who made helpful suggestions to improve the quality of the presentation. The work for this article,
in particular by the first, second, third, ninth, and last author, has been conducted within the Research Campus MODAL funded by the German Federal Ministry of Education and Research (BMBF grant numbers
05M14ZAM, 05M20ZBM). All responsibility is assumed by the authors. The work of Hans D. Mittelmann was supported in part by Air Force Office of Scientific Research under grants FA9550-15-1-0351 andFA9550-19-1-0070.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence,and indicate if changes were made. The images or other third party material in this article are includedin the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. Ifmaterial is not included in the article’s Creative Commons licence and your intended use is not permitted
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References
1. Achterberg, T.: Constraint integer programming. Ph.D. thesis, Technische Universität Berlin (2007)
2. Achterberg, T., Bixby, R.E., Gu, Z., Rothberg, E., Weninger, D.: Presolve reductions in mixed integer
programming. INFORMS J. Comput. (2019). https://doi.org/10.1287/ijoc.2018.0857
3. Achterberg, T., Koch, T., Martin, A.: MIPLIB 2003. Oper. Res. Lett. 34(4), 361–372 (2006). https://doi.org/10.1016/j.orl.2005.07.009
9 Compare the downloadable files open-v1.test and open-v13.test.

123

MIPLIB 2017: data-driven compilation...
4. Beasley, J.E.: OR-Library: distributing test problems by electronic mail. J. Oper. Res. Soc. 41(11),1069–1072 (1990)
5. Beyer, K., Goldstein, J., Ramakrishnan, R., Shaft, U.: When is “nearest neighbor” meaningful? In:
C. Beeri, P. Buneman (eds.) Database Theory–ICDT’99, pp. 217–235. Springer, Berlin (1999). https://doi.org/10.1007/3-540-49257-7_15
6. Bischl, B., Kerschke, P., Kotthoff, L., Lindauer, M., Malitsky, Y., Fréchette, A., Hoos, H., Hutter, F.,
Leyton-Brown, K., Tierney, K., et al.: ASlib: a benchmark library for algorithm selection. Artif. Intell.
237, 41–58 (2016)
7. Bixby, R.E.: Solving real-world linear programs: a decade and more of progress. Oper. Res. 50(1),
3–15 (2002). https://doi.org/10.1287/opre.50.1.3.17780
8. Bixby, R.E., Boyd, E.A., Indovina, R.R.: MIPLIB: a test set of mixed integer programming problems.
SIAM News 25, 16 (1992)
9. Bixby, R.E., Ceria, S., McZeal, C., Savelsbergh, M.: An updated mixed integer programming library:
MIPLIB 3.0. Optima 58, 12–15 (1998)
10. Browne, S., Dongarra, J., Grosse, E., Rowan, T.: The Netlib mathematical software repository. D-lib
Magazine 1(9), (1995)
11. Bussieck, M.R., Drud, A.S., Meeraus, A.: MINLPLib—a collection of test models for mixed-integer
nonlinear programming. INFORMS J. Comput. 15(1), 114–119 (2003)
12. COIN-OR branch-and-cut MIP solver (2019). https://github.com/coin-or/Cbc
13. IBM ILOG CPLEX Optimizer (2019). http://www-01.ibm.com/software/integration/optimization/
cplex-optimizer/
14. Crowder, H.P., Dembo, R.S., Mulvey, J.M.: Reporting computational experiments in mathematicalprogramming. Math. Program. 15, 316–329 (1978)
15. Fourer, R.M., Gay, D.W., Kernighan, B.: AMPL: A Modeling Language for Mathematical Programming. Duxbury Press, California (2002)
16. Furini, F., Traversi, E., Belotti, P., Frangioni, A., Gleixner, A., Gould, N., Liberti, L., Lodi, A., Misener,
R., Mittelmann, H., et al.: QPLIB: a library of quadratic programming instances. Math. Program.
Comput. 11(2), 237–265 (2019)
17. Galati, M., Ralphs, T., Wang, J.: Computational Experience with Generic Decomposition using the
DIP Framework. In: Proceedings of RAMP 2012 (2012). http://coral.ie.lehigh.edu/~ted/files/papers/
RAMP12.pdf
18. Galati, M., Ralphs, T., Wang, J.: DIP/DipPy version 0.92 (2019). https://doi.org/10.5281/zenodo.
2656800. http://github.com/coin-or/Dip
19. Galati, M.V.: Decomposition in Integer Programming. Phd, Lehigh University (2009). http://coral.ie.
lehigh.edu/~ted/files/papers/MatthewGalatiDissertation09.pdf
20. Gamrath, G., Koch, T., Martin, A., Miltenberger, M., Weninger, D.: Progress in presolving for mixed
integer programming. Math. Program. Comput. (2015). https://doi.org/10.1007/s12532-015-0083-5
21. Gamrath, G., Lübbecke, M.: Experiments with a generic Dantzig–Wolfe decomposition for integer
programs. In: P. Festa (ed.) Experimental Algorithms, Lect. Notes in Comp. Sci., vol. 6049, pp. 239–
252. Springer, Berlin (2010). https://doi.org/10.1007/978-3-642-13193-6_21
22. Georges, A., Gleixner, A., Gojic, G., Gottwald, R.L., Haley, D., Hendel, G., Matejczyk, B.: Featurebased algorithm selection for mixed integer programming. ZIB-Report 18-17, Zuse Institute Berlin(2018)
23. Gleixner, A., Bastubbe, M., Eifler, L., Gally, T., Gamrath, G., Gottwald, R.L., Hendel, G., Hojny,
C., Koch, T., Lübbecke, M.E., Maher, S.J., Miltenberger, M., Müller, B., Pfetsch, M.E., Puchert,
C., Rehfeldt, D., Schlösser, F., Schubert, C., Serrano, F., Shinano, Y., Viernickel, J.M., Walter, M.,
Wegscheider, F., Witt, J.T., Witzig, J.: The SCIP Optimization Suite 6.0. ZIB-Report 18–26. Zuse
Institute, Berlin (2018)
24. Granlund, T., the GMP development team: GNU MP: The GNU Multiple Precision Arithmetic Library,
6.1.2 edn (2016). http://gmplib.org/
25. GUROBI Optimizer (2019). http://www.gurobi.com/products/gurobi-optimizer/gurobi-overview
26. Hartigan, J.A., Wong, M.A.: Algorithm AS 136: A K-Means clustering algorithm. Appl. Stat. 28(1),100–108 (1979). https://doi.org/10.2307/2346830
27. Hendel, G.: IPET interactive performance evaluation tools. https://github.com/GregorCH/ipet
28. Hoffman, A., Mannos, M., Sokolowsky, D., Wiegmann, N.: Computational experience in solving linear
programs. J. Soc. Ind. Appl. Math. 1(1), 17–33 (1953)
29. Hooker, J.: Needed: an empirical science of algorithms. Oper. Res. 42(2), 210–212 (1993)

123

A. Gleixner, G. Hendel, G. Gamrath et al.
30. IBM Corporation, White Plains, N.Y.: Mathematical Programming System/360 Version 2, Linear and
Separable Programming—User’s Manual, 3 edn. (1969). Publication H20–0476–2
31. Jackson, R.H.F., Boggs, P.T., Nash, S.G., Powell, S.: Guidelines for reporting results of computational
experiments. Report of the ad hoc committee. Math. Program. 49, 413–425 (1991)
32. Koch, T., Achterberg, T., Andersen, E., Bastert, O., Berthold, T., Bixby, R.E., Danna, E., Gamrath, G.,
Gleixner, A.M., Heinz, S., Lodi, A., Mittelmann, H., Ralphs, T., Salvagnin, D., Steffy, D.E., Wolter,
K.: MIPLIB 2010. Math. Program. Comput. 3, 103–163 (2011). https://doi.org/10.1007/s12532-0110025-9
33. Laundy, R., Perregaard, M., Tavares, G., Tipi, H., Vazacopoulos, A.: Solving hard mixed integer
programming problems with Xpress-MP: A MIPLIB 2003 case study. INFORMS J. Comput. 21,
304–319 (2009)
34. Linderoth, J.T., Ralphs, T.K.: Noncommercial software for mixed-integer linear programming. Integ.
Program. Theory Pract. 3(253–303), 144–189 (2005)
35. Lodi, A.: MIP computation. In: Jünger, M., Liebling, T., Naddef, D., Nemhauser, G., Pulleyblank, W.,
Reinelt, G., Rinaldi, G., Wolsey, L. (eds.) 50 Years of Integer Programming 1958–2008, pp. 619–645.
Springer, Berlin (2009)
36. MathWorks: MATLAB Optimization Toolbox (2019). https://de.mathworks.com/products/
optimization.html
37. McGeoch, C.C.: Toward an experimental method for algorithm simulation. INFORMS J. Comput.
8(1), 1–15 (1996)
38. McGeoch, C.C.: Experimental analysis of algorithms. Notices AMS 48(3), 304–311 (2001)
39. MIPLIB 2.0. http://miplib2010.zib.de/miplib2/miplib2.html
40. Mittelmann, H.: Benchmarks for optimization software. http://plato.asu.edu/bench.html
41. MOSEK (2019). https://www.mosek.com/
42. Nazareth, J.L.: Computer Solutions of Linear Programs. Monographs on Numerical Analysis. Oxford
University Press, Oxford (1987)
43. Rice, J.R.: The algorithm selection problem. In: Advances in computers, vol. 15, pp. 65–118. Elsevier(1976)
44. SAS: SAS/OR (2019). https://www.sas.com/en_us/software/or.html
45. Solving Constraint Integer Programs (2019). http://scip.zib.de/
46. Shinano, Y.: The Ubiquity Generator framework: 7 years of progress in parallelizing branch-and-bound.
In: N. Kliewer, J.F. Ehmke, R. Borndörfer (eds.) Operations Research Proceedings 2017, pp. 143–149.
Springer International Publishing, Cham (2018). https://doi.org/10.1007/978-3-319-89920-6_20
47. Shinano, Y., Achterberg, T., Berthold, T., Heinz, S., Koch, T., Winkler, M.: Solving Open MIP Instanceswith ParaSCIP on Supercomputers using up to 80,000 Cores. In: Proceedings of the 30th IEEE International Parallel & Distributed Processing Symposium (2016). https://doi.org/10.1109/IPDPS.2016.
56
48. Smith-Miles, K., Baatar, D., Wreford, B., Lewis, R.: Towards objective measures of algorithm performance across instance space. Comput. Oper. Res. 45, 12–24 (2014)49. Smith-Miles, K., Bowly, S.: Generating new test instances by evolving in instance space. Comput.
Oper. Res. 63, 102–113 (2015)
50. van der Maaten, L., Hinton, G.: Visualizing high-dimensional data using t-SNE. J. Mach. Learn. Res.
9, 2579–2605 (2008)
51. Xpress: FICO Xpress-Optimizer (2019). http://www.fico.com/en/Products/DMTools/xpressoverview/Pages/Xpress-Optimizer.aspx
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps
and institutional affiliations.

123

MIPLIB 2017: data-driven compilation...

Aﬃliations
Ambros Gleixner1 · Gregor Hendel1 · Gerald Gamrath1 ·Tobias Achterberg2 · Michael Bastubbe3 · Timo Berthold4 ·Philipp Christophel5 · Kati Jarck6 · Thorsten Koch1 · Jeﬀ Linderoth7
Marco Lübbecke3 · Hans D. Mittelmann8 · Derya Ozyurt9 ·Ted K. Ralphs10 · Domenico Salvagnin11 · Yuji Shinano1·

Gregor Hendel
hendel@zib.de
Gerald Gamrath
gamrath@zib.de
Tobias Achterberg
achterberg@gurobi.com
Michael Bastubbe
bastubbe@or.rwth-aachen.de
Timo Berthold
timoberthold@fico.com
Philipp Christophel
Philipp.Christophel@sas.com
Kati Jarck
kati.jarck@mosek.com
Thorsten Koch
koch@zib.de
Jeff Linderothlinderoth@wisc.edu
Marco Lübbecke
luebbecke@or.rwth-aachen.de
Hans D. Mittelmannmittelmann@asu.edu
Derya Ozyurt
dozyurt@mathworks.com
Ted K. Ralphs
ted@lehigh.edu
Domenico Salvagnin
domenico.salvagnin@unipd.it
Yuji Shinanoshinano@zib.de
1

Zuse Institute Berlin, Takustr. 7, 14195 Berlin, Germany

2

Gurobi GmbH, Ulmenstr. 37-39, 60325 Frankfurt am Main, Germany

3

Lehrstuhl für Operations Research, RWTH Aachen University, Kackertstr. 7, 52072 Aachen,
Germany

4

Fair Isaac Germany GmbH, Stubenwald-Allee 19, 64625 Bensheim, Germany

5

Operations Research R&D, Advanced Analytics Division, SAS Institute Inc., Cary, NC, USA

123

A. Gleixner, G. Hendel, G. Gamrath et al.
6

MOSEK ApS, Fruebjergvej 3, Symbion Science Park, 2100 Copenhagen, Denmark

7

Department of Industrial and Systems Engineering & Wisconsin Institute for Discovery,
University of Wisconsin-Madison, Madison, WI, USA

8

School of Mathematical and Statistical Sciences, Arizona State University, Tempe, AZ
85287-1804, USA

9

The MathWorks, Inc., 1 Apple Hill Drive, Natick, MA 01760-2098, USA

10

Department of Industrial and Systems Engineering, Lehigh University, Bethlehem, PA 18015,
USA

11

Department of Information Engineering, University of Padova, Via Gradenigo 6/b,
35131 Padova, Italy123