1 Introduction
We live in a world full of tradeoffs. In almost every problem situation we encounter, it is difficult to define the one and only goal to aim for, especially whenever more than one decision maker or stakeholder is involved. Thus, many if not all practical problems involve several and often conflicting objectives. Prominent examples are profitability or cost versus environmental concerns (Ramos et al. 2014) or customer satisfaction (Braekers et al. 2016). Whenever it is not possible to aggregate these conflicting objectives in a meaningful way, it is recommendable to produce the set of tradeoff solutions in order to elucidate their tradeoff relationship.
In this paper, we propose a branchandbound framework which produces the optimal set of tradeoff solutions for optimization problems with two objectives which can be modeled as integer linear programs. The main contributions are a new problemindependent branching rule for biobjective optimization, building up on the Pareto branching idea of
Stidsen et al. (2014), and several improvements exploiting the integrality of objective function values. We compare our framework to existing general purpose stateoftheart frameworks, such as the constraint framework and the balanced box method, and to results obtained by means of another branchandbound algorithm developed by Belotti et al. (2013). Furthermore, we apply it to a biobjective orienteering problem for which lower bound sets are produced using column generation. Although column generation has been used in the context of biobjective optimization, e.g. in the context of robust airline crew scheduling (Tam et al. 2011), vehicle routing (Sarpong et al. 2013), and biobjective linear programming (Raith et al. 2012), this is, to the best of our knowledge, the first time column generation is used in a biobjective branchandbound algorithm, resulting in a biobjective branchandprice algorithm.The paper is organized as follows. In Sections 2 and 3, we give preliminaries and we review related work. In Section 4, we present the key ingredients to our branchandbound framework and in Section 5 we discuss the proposed enhancements. Experimental results and how we embed column generation into our branchandbound framework are presented in Section 6. Section 7 concludes the paper.
2 Preliminaries
The aim of this paper is to solve biobjective integer linear programs of the following form:
s.t.  (1)  
where
is the vector of objective functions,
is the constraint matrix of size , the right hand side vector of size , and the decision variables which are restricted to integer values. An important subclass are binary (combinatorial) problems.In formulation (1) we strive for identifying the set of Pareto optimal solutions according to the given objectives. A feasible solution is Pareto optimal if there does not exist any feasible solution which dominates it, i.e. which is better in one objective and not worse in any other objective. The image of a Pareto optimal solution in criterion space is referred to as nondominated point. Each nondominated point may represent more than one Pareto optimal solution. We aim at generating all nondominated points, i.e. at least one Pareto optimal or efficient solution per nondominated point.
Depending on the structure of the objective functions and the feasible domains of the variables, the set of all nondominated points may take different shapes. In the case where all variables assume continuous values, it corresponds to the boundary of the convex hull of the image of all feasible solutions in criterion space. It can be characterized by the extreme points of the boundary. These extreme points can be computed efficiently, e.g. by the weighted sum method of Aneja and Nair (1979). It combines the two objectives into a single one, resulting in , with and giving the weights of objective one and two, respectively. These weights are changed systematically so as to obtain the entire set of extreme points. In the case where all variables may only take integer values, which we look at, the set of nondominated points is composed of points which may or may not lie on the boundary of the convex hull in criterion space. Solutions whose images in criterion space correspond to extreme points are called extreme (supported) efficient solutions, those that correspond to points that lie on the boundary are called supported efficient solutions and those that correspond to points that lie in the interior of the convex hull but are still nondominated are called nonsupported efficient solutions. In the integer case, the set of extreme points can also be computed by means of the weighted sum method.
Additional important notions in multiobjective optimization are ideal point and nadir point. The ideal point is constructed by combining the best possible values for each objective function. It dominates all nondominated points. The nadir point is its opposite, a point that is dominated by all nondominated points. It is constructed by combining the worst possible values for each objective function of all efficient solutions. Both notions have been used in the context of branchandbound algorithms for multiobjective optimization. For a more detailed introduction to multicriteria decision making, we refer to Ehrgott (2005).
3 Related work
Exact approaches in multiobjective (mixed) integer programming can be divided into two classes: those that work in the space of objective function values (referred to as criterion space search methods, e.g. by Boland et al. (2015a)) and those that work in the space of decision variables (generalizations of branchandbound algorithms).
Criterion space search methods solve a succession of singleobjective problems in order to compute the set of Pareto optimal solutions. Therefore, they exploit the power of stateoftheart mixed integer programming solvers. This appears to be one of their main advantages in comparison to generalizations of branchandbound algorithms (Boland et al. 2015b). However, many combinatorial singleobjective problems cannot be solved efficiently by commercial solvers, e.g. the most efficient exact algorithms in the field of vehicle routing rely on column generation based techniques (see e.g. Baldacci et al. 2012). Thus, especially in this context but also in general, it is not clear whether a criterion space search method or a biobjective branchandbound algorithm is more efficient.
One of the most popular criterion space search methods is the constraint method, first introduced by Haimes et al. (1971). It consists in iteratively solving singleobjective versions of a biobjective problem. In every step, the first objective is optimized, but a constraint is updated in order to improve the quality of the solution with regards to the second objective. Thus all nondominated points are enumerated. The constraint method is generic and simple to implement and it is among the best performing criterion space search algorithms when applied, e.g. to the biobjective prizecollecting Steiner tree problem (Leitner et al. 2014).
Recently, the balanced box method for biobjective integer programs has been introduced by Boland et al. (2015a). Optimal solutions are computed for each objective and they define a rectangle. This rectangle is then split in half and in each half, the objective which has not been optimized yet in this half of the rectangle is then optimized. At this stage there are two nonoverlapping rectangles. The same procedure is repeated recursively on these rectangles.
Generalizations of branchandbound to multiple objectives for mixed 01 integer programs are provided by Mavrotas and Diakoulaki (1998). Their bounding procedure considers an ideal point for the upper bound (they work on a maximization problem), consisting of the best possible value for each objective at this node, and keeps branching until this ideal point is dominated by the lower bound set. Vincent et al. (2013) improve the algorithm by Mavrotas and Diakoulaki (1998), most notably by comparing bound sets instead of ideal points in the bounding procedure. Masin and Bukchin (2008) propose a branchandbound method that uses a surrogate objective function returning a single numerical value. This value can be treated like a lower bound in the context of single objective minimization and conveys information on whether the node can be fathomed or not. They illustrate their method using a threeobjective scheduling problem example but no computational study is provided. Sourd and Spanjaard (2008)
use separating hyperplanes between upper and lower bound sets in order to discard nodes in a general branchandbound framework for integer programs. The concept is in fact similar to the lower bound sets defined by
Ehrgott and Gandibleux (2006). Jozefowiez et al. (2012) introduce a branchandcut algorithm for integer programs in which discrete sets are used for lower bounds, so nodes can be pruned if every point in the lower bound set is dominated by a point in the upper bound set. Partial pruning is used in order to discard parts of the solutions represented by a given node of the search tree, thus speeding up the whole tree search. Belotti et al. (2013) are the first to introduce a generalpurpose branchandbound algorithm for biobjective mixed integer linear programming, where the continuous variables may appear in both objective functions. They build up on the previous work by Visée et al. (1998), Mavrotas and Diakoulaki (1998), Sourd and Spanjaard (2008) and Vincent et al. (2013) and like them, they use a binary branching scheme. Improved fathoming rules are introduced in order to discard more nodes during tree search. Adelgren et al. (2014) propose an efficient data structure to store and update upper bound sets in the context of mixed integer programming and illustrate its efficiency using the algorithms of Belotti et al. (2013) and of Adelgren and Gupte (2016).Ideas most related to ours are presented by Stidsen et al. (2014). They provide a branchandbound algorithm to deal with a certain class of biobjective mixed integer linear programs. They propose two binary branching schemes that exploit partial dominance of already visited integer solutions and the current upper bound set, respectively. Our approach relies on bound sets to represent the lower bound, in contrast to Stidsen et al. (2014) who solve one scalarized single objective problem in each node of their branchandbound tree. In addition, in our scheme, the number of child nodes depends on the number of currently nondominated portions which are obtained after applying what we call filtering of the lower bound set by the current upper bound set. Ideas similar to the Pareto branching concept of Stidsen et al. (2014) and its extension to bound sets presented in this paper, which we name objective space branching, are also being explored by Adelgren and Gupte (2016) in the context of a branchandbound method for multiobjective mixed integer programming and by Gadegaard et al. (2016a) in the context of a biobjective branchandcut algorithm.
Furthermore, we also propose several improvements exploiting the integrality of objective functions. These improvements, which we call segment tightening, lower bound lifting, and integer dominance allow to disregard even larger regions of the objective space. We illustrate their impact on the overall performance of the algorithm in the experimental results section.
4 A branchandbound framework for biobjective integer optimization
Branchandbound is a general purpose tree search method to solve (mixed) integer linear programs. As its name suggests, it has two main ingredients: branching and bounding. Branching refers to the way the search space is partitioned, while bounding refers to the process of using valid bounds to discard subsets of the search space without compromising optimality. In the following, we first describe the general framework of our biobjective branchandbound (BIOBAB). We then describe its different ingredients in further detail.
4.1 General idea
Our BIOBAB is similar to a traditional singleobjective branchandbound algorithm, except for the fact that we compare bound sets, instead of single numerical values.
In the beginning the upper bound (UB) set is empty. It is updated every time an integer solution is obtained. One of three things can happen:

if is dominated by at least one solution from the UB set, then is discarded,

if some solutions from the UB set are dominated by , then these solutions are discarded and is added to the UB set,

if neither dominates nor is dominated by any solution from the UB set, is directly added to the UB set.
These conditions can be tested in linear time. If the UB set is sorted according to the values of one objective, it is possible to check dominance in logarithmic time using binary search.
In each branchandbound node, lower bound (LB) sets are computed by means of an algorithm similar to the weighted sum method of Aneja and Nair (1979), described in detail in Section 5.2, with either the linear relaxation or the original integer program (IP), since both yield valid bound sets (Ehrgott and Gandibleux 2006).
A branchandbound node contains information on branching decisions as usual, but also a bound for each objective. This is used to partition the objective space. The objective bounds at the root node can be known valid bounds, or infinite values.
Once the LB set at a given node has been calculated, we apply what we call filtering (see Section 4.3): the current UB set and the current LB set are compared and a set of nondominated portions of the current LB set are returned. In the case where this set is empty, the node can be fathomed. Otherwise branching may be required; the filtered LB set is then used as input for the branching procedure (see Section 4.4). If the LB set corresponds to a leaf (i.e. the node corresponds to a unique integer feasible solution), or if it is completely dominated following filtering, then the node can be fathomed and branching is not necessary. Otherwise a number of child nodes are generated, which are added to the list of branchandbound nodes to be processed. Once all nodes have been explored, the search is over and the complete Pareto set for the original problem is contained in the UB set.
We now present the concept of lower bound segment which is used in the bounding, filtering and branching procedures and which is also key to the proposed improvements discussed in Section 5.
4.2 Lower bound segments
In the following we assume that a LB set has been produced. As mentioned above, it is described by the extreme points of the convex hull (of the image of the feasible set in objective space) for the linear relaxation of the problem (or the original IP).
To any two consecutive extreme points and we can associate the subset of objective space that is covered or dominated and that represents the possibility to find feasible solutions to the integer problem. To define this covered space, we associate to each pair of extreme points a topright corner (or local nadir) point of which coordinates are valid singleobjective upper bounds, one per objective, on the space dominated by the two extreme points and the line connecting them.
Let be the set of all points in the objective space which are associated to feasible solutions to the biobjective integer problem at hand. We now formally introduce the concept of lower bound segment that we use to represent LB sets.
Definition 1
Three points , and , such that , define a lower bound segment iff , where is the slope of the line defined by the two points and is the yintercept of this same line; is a valid local nadir point, such that and .
Thus, any point is on or above the line defined by and . Note that we do not include the case where into our definition. For convenience, we extend the definition of segments to allow equal points:
Definition 2
Two points and define a lower bound segment iff there does not exist a feasible point that dominates ; is a valid local nadir point, such that and
Figure 1 depicts such a lower bound segment (shaded).
In the following we use objectoriented notation to refer to the attributes defining a given segment: if extreme points , and local nadir point define a segment , then, , and . Additionally, is the slope of the line going through and , while is its intercept. Graphically, a segment is a convex polygon defined by the points , , , , .
In the most basic case the coordinates of the local nadir point can be arbitrarily large values. If the two extreme points of the set of nondominated points are known, we can deduce a tighter valid local nadir point from them. Similarly, we can associate to any LB set a valid local nadir point by considering, for each objective, the maximum value among all local nadir points.
Here we note that we can partition the objective space, using values for any of the two objectives, or any linear combination of both with positive weights (as is done in Stidsen et al. 2014). In our BIOBAB, we partition the objective space using different intervals for the first objective to split the LB set into vertical stripes.
In practice, this allows us to improve the initial local nadir point for segments of a connected sequence of segments: for any two segments connected by point , we can cut from the space covered by the left segment all points such that , since these points are also covered by the right segment. By doing so, we partition the objective space covered by the LB set into different LB segments with associated local nadir points.
UB sets can also be used to provide better local nadir points. For instance, following our previous example from Figure 1, consider point , being the image of a feasible solution in objective space, such that and , then defines a valid bound on the first objective of the nadir point, since any point such that and is dominated by . This example is described graphically in Figure 2: for the depicted segment, using UB point , the initial local nadir point can be improved to .
The fact that LB sets can be reduced using UB sets is the basis for a biobjective branching rule that is presented in Section 4.4.
4.3 Bound set filtering and node fathoming
One major difficulty in previous biobjective branchandbound approaches lies in the evaluation of the dominance of a given LB set by a given UB set. Multiple fathoming rules have been developed over the years (see Belotti et al. 2016, for the current state of the art). We now introduce the fathoming rule used in BIOBAB.
As seen above, it is possible to split a LB set into LB segments (with associated local nadir points). If we establish that each segment of a current LB set is dominated by the UB set, then we also establish that the whole LB set is dominated, therefore the current node of the branchandbound tree can be fathomed. In order to determine if a LB segment is dominated by an UB set, we simply subtract the space dominated by the UB set from the LB segment: if the remaining space is empty then the LB segment is dominated, otherwise this segment may be tightened (as shown for example in Figure 2).
Algorithm 1 details how to subtract the space covered by an UB point from a LB segment . It shares some similarities with the dominance rules used by Vincent et al. (2013), but considers more cases. Notably, it updates the local nadir point even in cases where does not dominate any point on the line connecting and . We use the objectoriented notation defined above. A pair represents a point and is used to create a segment out of points and . Lines 27 deal with the case where does not dominate any point on the line connecting and . In that case, may still be above or to the right of , in which cases there is potential for improving . If is to the right of and below (Lines 57), we are in the situation depicted in Figure 2 and can be reduced to . Similarly, if is above and to the left of (Lines 34), can then be reduced to . Lines 815 deal with all the cases when at least one part of the line connecting and is dominated by . There are four different cases: (i) does not dominate (Lines 911), (ii) does not dominate (Lines 1214), (iii) dominates neither nor , resulting in two new segments (both clauses in Lines 911 and 1214 are matched), and (iv) dominates both and . In case (iv) no clause is matched in the algorithm, therefore an empty set is returned.
Now in order to filter a whole LB set using an UB set, we filter each segment in the LB set with each point in the UB set. If the remaining LB set is empty then the node can be fathomed, otherwise branching is required.
4.4 Branching procedure
When dealing with biobjective bound sets, a major challenge is that the LB set at a given node of the branchandbound tree may be partially dominated but not completely, therefore the node cannot be fathomed. This is illustrated in Figure 3.
In order to disregard those portions of the objective space which are dominated by the UB set, we use the following branching rule which we name objective space branching: each contiguous subset of the filtered LB set gives rise to a new branch and the objective space between those parts is simply ignored. To generate the branch for a given subset, we consider its local nadir point and add two cutting planes so that both objectives are better than the value at this local nadir point: and . In practice, this is achieved by setting the objective bounds of the corresponding child nodes using .
Branching on objective space has no impact on decision space, i.e. computing lower bounds for these child nodes will generate the same LB set, albeit in several pieces, one per branch in objective space. Therefore, objective space branching (OSB) never happens alone; branching on decision space is systematically also performed at the same time. The end result is that each node in the search tree contains at least one new branching decision on decision space, and optionally a branching decision on objective space.
In general, decisionspace branching is problem specific. However, we can still make the general remark that there is no reason to systematically use the same variable for decisionspace branching in every OSB child node: each OSB child has its own separate LB set, and there is no reason to branch on the same decision variables when considering different disjoint LB sets. For this reason, OSB is always applied first in order to reduce the search space; for each OSB child and associated LB set, decisionspace branching is conducted independently. The end result might still be that each OSB child branches on the same variable, depending on the branching rule.
Algorithm 2 gives the function of our branchandbound framework. It takes a node of the branchandbound tree and an LB set as input and returns a set of nodes. A node is represented as a pair where represents the objective bounds (local nadir point) for this node and is the set of branching decisions for this node.
In the example of Figure 3, function considers four disjoint regions, each of them generating their own set of branching decisions.
We note here that even if an IP is used to compute the LB set, any given variable may take different values over a same LB set; therefore branching on decision space is still necessary.
5 Improvements based on the integrality of objective functions
When solving integer programs, in many cases, objective values of feasible integer solutions only take integer values. For pure integer programs, it is in practice almost always possible to use integer coefficients. The reasons include the fact that LP solvers have precision limitations, and that timeefficient floatingpoint numbers representations also have precision limitations. So rounding floatingpoint numbers is almost always inevitable, and if numbers are rounded to decimals then they may as well be multiplied by and considered integers. We note that this property is exploited by other methods as well. For instance, the constraint framework uses a known value which in our case would be 1. The balanced box method also relies on similar values. In some cases, it is even possible to use values higher than 1, as long as these values are valid: for instance if every coefficient for a given objective function is integer, then the greatest common divisor of these coefficients can be used as a valid value. For the sake of simplicity and without loss of generality, we consider in the following that this valid value is 1. We now explore possibilities to enhance our biobjective branchandbound by exploiting the fact that objective values of feasible solutions are always integer.
Any given LB segment covers a part of the objective space, but also includes subsets not containing any point with integer coordinates. Such subsets can be disregarded during the search. This is illustrated in Figure 4, where dashed lines represent integer values for each objective.
This general idea can be exploited in several ways in order to speed up the search.
5.1 Segment tightening
First, any LB segment which does not cover any integer vector can be discarded. This can be tested in .
We first investigate the case where , then and are defined:
Proposition 1
Segment with covers integer vectors iff .
Proof 1
A segment is a polygon defined by the points , , , , . The coordinates of any integer point within this polygon have to satisfy and and . The integer point closest to the local nadir point which has the potential to lie within the polygon has coordinates . Any other integer point contained within the region of the segment must have coordinates and . Therefore, if is contained within , covers at least one integer point.
The special case of segments defined by two points only (see Definition 2) is even simpler, as the space covered by such a segment is a rectangle. Consider the segment defined by point and local nadir point , then this segment contains integer points iff .
This is tested right after filtering: only the segments in the LB set that cover an integer point are kept. This enhancement is called segment tightening.
5.2 Lower bound lifting
Second, we can use the fact that feasible integer solutions have integer objective values in order to speed up our LB set computation procedure. We call this procedure (it is given in Algorithm 3).
As in the original algorithm of Aneja and Nair (1979), procedure successively computes supported points, starting with the two points at the two ends of the efficient frontier, using lexicographic minimum objective functions and . This is similar to what is done by Boland et al. (2015a) in the context of the balanced box method.
For any two given points, solving a certain weightedsum problem determines whether there exists another supported solution “between” these two. However, it can happen that these two points already describe a valid LB segment (given an appropriate nadir point), regardless of any other extreme supported point between them. If there is no point with integer coordinates in the triangle defined by , and , then and already define a valid LB segment and there is no need to investigate these two points further. We note that point can also be disregarded because it dominates both and , although it is already established that they are both nondominated. Moreover, if the line connecting and is not a part of the convex hull boundary then it defines a tighter lower bound than the one induced by the actual convex hull boundary between and . This situation is illustrated in Figure 5, where and define a valid lower bound segment (given an appropriate local nadir point).
As previously, testing for the existence of an integer vector in the given triangle can be achieved in and sometimes allows to skip some stages of the algorithm while providing a tighter LB set. In the notation of Algorithm 3, this means that after retrieving from (line 10, Algorithm 3), we check if the triangle derived from does not contain any integer point.
Proposition 2
If the point is above the line connecting and , i.e. , then the triangle derived from does not contain any integer point.
Proof 2
Since and are extreme points of the boundary of the convex hull, an additional LB point may only be found in the right triangle defined by the points , , , with and . Let us call the local ideal point. Since and , must either lie within the triangle or it is above the line connecting and . Since any other integer point must have coordinates satisfying and , if lies above the line connecting and , no other integer point within this triangle exists.
In this case, the segment derived from is appended to in Algorithm 3 and lines 13–21 are skipped. Otherwise a weighted sum problem () is solved, where and are determined by the algorithm. In the case where the resulting point lies below the line connecting and , and are stored in for further processing, as in the algorithm by Aneja and Nair (1979). This enhancement is called lower bound lifting.
5.3 Integer dominance
Third, it is also possible to exploit the integrality of objective values for feasible integer solutions when branching on objective space. If we are in the situation of branching on objective space, it is typically because part of the LB set has been dominated by the UB set. In that case, we already know that the given local nadir point of any LB segment is in fact also dominated by the UB set (this is actually how this local nadir point is calculated, see Algorithm 1). So this local nadir point can be tightened.
Figure 6 illustrates this principle: after filtering the LB segment with points from the upper bound set, and (Figures 6 and 6), we obtain a reduced segment . The new local nadir point is actually dominated by both and (Figure 6). Moreover, each point in the area which is in a darker shade (Figure 6) is either (i) dominated by and/or or (ii) not integer. Therefore we can disregard the whole area. This can be achieved by subtracting any value such that from both coordinates of the local nadir point after filtering, thus producing an improved segment . In order to implement this idea in our algorithm, we simply subtract from both coordinates of each upper bound point used in Algorithm 1. This enhancement is called integer dominance. It is in fact similar to using an value of 1 in the constraint framework, or in the balanced box method.
6 Computational study
We investigate the efficiency of the proposed algorithm and enhancements through experimentation.
A generic version of BIOBAB is developed in Python 2.7 and Gurobi 6.5.0 is used to solve LPs and MIPs. Algorithms are run on an Intel Xeon E52650v2 CPU at 2.6 GHz, with a 4 GB memory limit and a twohour CPU time limit. Except when explicitly stated otherwise, BIOBAB uses the linear relaxation of the original IP to compute lower bound sets. In terms of tree exploration, BIOBAB always performs a breadthfirst search. Like in Boland et al. (2015a), we use performance profiles to compare the performances of different algorithms (Dolan and Moré (2002)). The performance of an algorithm for a certain instance is the ratio of the CPU time required by that algorithm for that instance to the best CPU time for that instance. The axis represents performance, while the axis represents the fraction of instances solved at a certain performance or better. If an algorithm does not terminate (i.e. find the Pareto set) for a certain instance, then it does not offer a performance for that instance. If no algorithm terminates for a given instance, then this instance is discarded.
In the following, we first describe the considered problems and benchmark instances for the first set of experiments, where we evaluate the different components of the algorithm and the impact of the integralitybased improvements. We then compare the results of our method with those of criterion space search algorithms and to an existing branchandbound method. Thereafter, we formally introduce the biobjective team orienteering problem with time windows (BITOPTW), whose linear relaxation is solved by means of column generation, we explain the different components tailored to this problem, and we report the obtained results.
6.1 Considered problems and benchmark instances
For our first set of experiments, we apply BIOBAB to two different uncapacitated biobjective facility location problems. The first uncapacitated facility location problem, which we call uncapacitated biobjective facility location problem (UBOFLP), considers the cost of opening facilities as first objective and the coverage of the population as second objective. Population from a certain node can only be covered by an open facility if the distance between the two is below a threshold. For the UBOFLP, we use the data that was also used to generate instances from the biobjective stochastic covering tour problem presented in Tricoire et al. (2012). The following simplifications are performed: (i) the transportation costs are ignored, (ii) the stochastic aspect is ignored and (iii) the distance aspect is simplified so that a node is either entirely covered or not covered at all by an open facility. In order to produce instances with various levels of difficulty, we concatenate the original instances: the UBOFLP instance consists of the data of the first covering tour problem instances. In total there are 36 UBOFLP instances, including 21 to 500 locations. The second uncapacitated facility location problem we consider is a simplification of the singlesource capacitated facility location problem, where the capacity constraints are relaxed. The first objective is to minimize the cost of opening facilities, while the second objective is to minimize assignment costs. We name this problem SSUFLP from now on. We use the 120 instances publicly available from Gadegaard et al. (2016b)
. For both UBOFLP and SSUFLP, decision space branching is performed in a generic way. First, the average value of each binary variable over the whole LB set (or rather subset, see Section
4.4) is computed. Then, the variable with average value closest to 1 is selected for binary branching on that subset. We provide mathematical formulations for UBOFLP and SSUFLP in Appendices B and C. For each of these problems, we compute the greatest common divisor for each objective and use it as valid value in the various methods. Since all costs are integer, is always at least 1.We provide information on the cardinality of efficient sets for UBOFLP and SSUFLP instances in Figure 7.
6.2 Impact of integralitybased improvements
In a first series of experiment, we assess the efficiency of the improvements described in Section 5, using the UBOFLP and SSUFLP data set. We run a version of BIOBAB with no improvement at all, a version with objective space branching only, a version with segment tightening only, a version with lower bound lifting only, and a version with all three improvements. Performance profiles are displayed in Figure 8.
On UBOFLP instances, full BIOBAB as well as segment tightening clearly outperform the other variants. Whether full BIOBAB outperforms the version with segment tightening only is unclear, because the whole chart is scaled to also include the performance profiles of poorly performing variants. For the same reason, it is also unclear whether lower bound lifting brings a significant improvement over bare BIOBAB. On SSUFLP instances, the differences are smaller but still clearly marked; full BIOBAB is still the best, but this time it appears that objectivespace branching is working better on its own than segment tightening. Again, it is unclear whether lower bound lifting brings any improvement.
We further investigate the individual contributions of the various BIOBAB components, by comparing full BIOBAB to full BIOBAB minus individual improvements. We first look at segment tightening, which appears to be crucial in certain cases, as shown in Figure 8. Figure 9 presents the performance of full BIOBAB with and without segment tightening.
Approximately 70% of UBOFLP instances cannot be solved within two hours when segment tightening is disabled, which makes it a crucial component. All SSUFLP Instances can be solved, but the profile of full BIOBAB is still consistently better than the profile without segment tightening. Full BIOBAB is up to 56% faster.
We now also look at the performance profile of full BIOBAB, versus full BIOBAB minus lower bound lifting, depicted in Figure 10.
The impact of lower bound lifting is less spectacular, but still significant. In rare cases not using lower bound lifting actually improves performance, but only very slightly (less than 5%, on some UBOFLP instances). Lower bound lifting speeds up the algorithm by up to 34% on UBOFLP instances and by up to 60% on SSUFLP instances. As mentioned earlier, lower bound lifting only involves constant time operations, so the overhead is minimal and will never cause an important decrease in performance, even when it is entirely useless. Therefore it is safe to use it in general.
Full BIOBAB is compared to the variant without OSB in Figure 11.
The impact of objective space branching is important: using it speeds up BIOBAB by a factor up to 14.7 for UBOFLP instances and up to 2.06 for SSUFLP instances. We note here that it is not surprising that objective space branching works well when combined with segment tightening, since segment tightening creates discontinuities in LB sets, thus allowing objective space branching to be applied.
These performance profiles clearly show that every algorithmic improvement brings added value, with rare exceptions for lower bound lifting. In order to provide additional insight on the use of objective space branching, we consider additional information related to full BIOBAB and full BIOBAB minus OSB: we analyze CPU time, number of nodes explored during tree search and number of linear programs solved during the whole search.
Full BIOBAB  No obj. space branching  
Instance  # LPs  # nodes  CPU (s)  # LPs  # nodes  CPU (s) 
021  23  1  0.30  23  1  0.28 
044  59  3  0.57  59  3  0.61 
055  73  3  0.89  73  3  0.85 
072  95  3  1.59  95  3  1.56 
090  135  7  3.04  197  7  4.18 
106  157  7  5.07  235  7  6.33 
120  173  7  6.44  267  7  8.88 
132  193  9  8.30  335  9  13.52 
134  197  9  8.58  339  9  15.16 
135  199  9  8.41  333  7  14.72 
162  263  11  16.78  401  9  26.39 
182  466  17  32.53  699  13  51.58 
203  478  17  41.57  697  11  66.53 
226  255  9  31.41  319  9  44.19 
254  259  9  41.72  335  9  53.46 
264  269  9  47.21  345  9  58.12 
266  271  9  49.48  347  9  63.44 
275  277  9  53.35  303  7  56.18 
294  285  9  64.19  361  9  84.88 
295  277  7  71.97  363  9  88.12 
296  289  9  74.64  365  9  94.16 
326  361  15  112.16  1459  29  467.91 
342  399  19  137.47  1125  15  397.68 
355  473  31  165.05  3509  61  1379.35 
370  455  25  195.29  983  19  386.76 
388  461  21  204.75  3159  39  1494.41 
399  531  35  236.88  1449  33  630.14 
410  595  45  297.69  1061  17  551.82 
428  645  51  342.32  7453  89  3885.49 
436  609  43  351.86  1585  31  919.27 
449  679  53  400.70  10586  183  5891.94 
458  713  57  432.92  5387  85  3241.36 
472  719  59  462.29  4681  113  2790.82 
482  821  75  591.75  5511  81  3456.99 
499  773  67  542.06  4209  97  2912.71 
500  759  65  583.60  9003  135  6927.28 
Table 1 shows the perinstance difference between the two versions of BIOBAB, on UBOFLP instances. We note that the total number of nodes explored does not vary that much between both methods, but the number of LPs solved does: the version without objective space branching solves more than ten times more LPs in some cases. We believe this validates the idea that branching on objective space avoids visiting several times regions that have already been established to be dominated by the UB set; apparently visiting these regions requires solving many LPs.
We now look closer at the number of branches generated by OSB. For that purpose we count, for each instance solved, how many branches are generated by OSB. We differentiate the number of branches generated at the root node and in the remainder of the branchandbound tree. We report average values for these two indicators, as well as the maximum tree depth at which OSB happened, in Table 2. Instances that are entirely solved without branching are disregarded.
Problem 




SSUFLP  23.84  2.07  12.24  
UBOFLP  5.37  0  0 
# root: number of OSB branches at the root node
# other: number of OSB branches at other nodes
max. depth: maximum depth at which OSB occurred
For both problem classes, OSB generates more branches at the root node, thus indicating that integer solutions are already found while solving the root node. For UBOFLP, OSB actually only happens at the root node; for SSUFLP, the average number of branches at nodes other than the root node is very close to 2, which suggests that it typically happens when one new efficient point is found. We also observe that a lower number of branches at a given node does not necessarily indicate that fewer integer solutions are found at this node, since segment tightening can also eliminate the objective space between two integer points, thus reducing the number of branches. But if there are OSB branches following the bounding procedure at a given node, then we know that at least integer solutions were found during this bounding procedure.
6.3 Comparison with criterion space search algorithms
We implemented a generic balanced box method as described in Boland et al. (2015a) (it is sketched in Algorithm 5 in the Appendix), including solution harvesting and other enhancements, within the same Python code base as BIOBAB. The frequency with which solution harvesting is called is controlled using parameter . Like in the original paper, we set following preliminary experiments. Also like in the original paper, we find that in order to produce lexicographic minima, solving explicitly two successive IPs is more efficient than solving one weightedsum IP. Previously found integer solutions that are valid for the problem at hand (i.e. within the same rectangle) are used to provide a cutoff value to the MIP solver.
Within the same code base we also develop generic versions of the constraint framework (see Algorithm 4 in the Appendix) as well as the bidirectional constraint framework described in Boland et al. (2015a). We note here that there are always two ways to run the singledirectional constraint algorithm, either starting with objective 1 or with objective 2. For a given problem class, one of these two ways is typically faster than the other, sometimes much faster. In these cases, it is likely that the bidirectional version will not be as fast since it has to solve IPs in the slower direction for 50% of these IPs. All constraint methods benefit from the applicable enhancements developed for the balanced box method.
We first compare BIOBAB and other criterion space search methods on the UBOFLP. We note here that BIOBAB solves the linear relaxation of the UBOFLP to compute LB sets, whereas criterion space search methods call the MIP solver directly. Performance profiles are given in Figure 12.
BIOBAB appears to perform better than criterion space search methods on the UBOFLP. It is however not overwhelming, as all methods always converge within a factor 3 of the best performance overall, i.e. within the same order of magnitude. Nonetheless, BIOBAB is definitely competitive on these instances. We also note that singledirectional constraint performs better than its bidirectional counterpart.
We now look at the performance profiles of the same methods, this time on SSUFLP instances. These profiles are represented in Figure 13.
The picture is now different: BIOBAB is the worst of all compared methods. It is however, still in the same order of magnitude (7 times slower than the best method in the worst case). In this comparison, all criterion space search methods benefit from the power of commercial MIP solvers, which are especially efficient at tree exploration and branching on decision variables. Even though we show above that objective space branching significantly improves the performance of BIOBAB, branching on decision variables is still performed in a somewhat naive way (fractional variable with value closest to 1, breadthfirst tree exploration). In the next experiment we do not relax integrality constraints when computing LB sets for the SSUFLP: IPs are solved within the bounding procedure of BIOBAB. An expectation is that LB sets will be of better quality, at the cost of extra CPU time. Better LB sets should reduce the size of the branchandbound tree of BIOBAB, thus delegating part of the task of tree search to the commercial MIP solver. Another expectation is that more integer solutions will be obtained earlier in the search, which in turn will allow OSB and integer dominance to speed up the search. We now compare this version of BIOBAB to the same criterion space search methods, in Figure 14.
The performance profiles clearly show that using an IP instead of its linear relaxation can be beneficial: BIOBAB is now performing much better, actually slightly better than all criterion space search methods.
Additionally, we display in Figure 15 the total number of linear programs (LPs) solved by each algorithm to generate the whole set of nondominated points. It is clear that the poor performance of BIOBAB using LPs is correlated with the unusually high amount of LPs it requires to solve, while BIOBAB using MIPs requires to solve a number of LPs comparable with that of other algorithms. This strengthens the idea that efficient branching on decision variables is the issue here.
6.4 Comparison with biobjective branchandbound
We now compare BIOBAB with the branchandbound algorithm for biobjective mixedinteger programs from Belotti et al. (2013). Their algorithm is for mixedinteger programs, but they also provide results on the biobjective set covering problem, which only uses binary variables. We run BIOBAB on the same instances and can therefore compare the CPU time required for solving these instances. We present results using the same format as in their article. Their computer is a 3.2 GHz workstation with 4 GB of RAM. We observe that BIOBAB is always faster, more than 10 times faster for the largest instances. Even if it is not clear how much has to attributed to differences in CPU, we can safely conclude that BIOBAB can compete with a stateoftheart branchandbound method on these instances.
Instance size ()  # instances  CPU (s) Belotti et al.  CPU (s) BIOBAB 
5 25  1  0.1  0.01 
10 50  1  2.1  0.73 
10 100  4  49.3  3.76 
Belotti et al.: 3.2 GHz workstation with 4 GB of RAM.
BIOBAB: 2.6 GHz Xeon CPU with a 4 GB RAM limit.
6.5 Application to cases which cannot benefit from the power of MIP solvers
A major strength of stateoftheart MIP solvers lies in their ability to conduct efficient branchandbound or branchandcut tree search, through involved tree exploration heuristics that are industrial secrets. However, it is not always possible to benefit from these advanced methods. Sometimes even the single objective problem version is too hard to solve with a MIP solver. We want to assess how BIOBAB fares in such cases. For that purpose, we use the BITOPTW as test case. The BITOPTW is formally described in Appendix
D. Such a model is challenging for current stateoftheart MIP solvers.6.5.1 Lower bound set: column generation
Since stateoftheart exact methods for singleobjective routing problems mostly rely on column generation based techniques (see, e.g. Baldacci et al. 2012), we also generate lower bound sets for the BITOPTW by means of column generation.
Let denote the total score or profit achieved by route , the total travel cost of route , and let indicate whether location is visited by route () or not (). Using binary variables equal to if route is selected from the set of all feasible routes, the BITOPTW can also be formulated as a pathbased model:
(2)  
(3) 
(4)  
(5)  
(6) 
Relaxing integrality requirements on the variables, we replace constraints (6) with:
(7) 
We also combine the two objective functions into a weighted sum ( and giving the respective nonnegative weights determined by our LB set computation procedure ):
(8) 
We thus obtain a single objective linear problem that can be solved by means of column generation, which allows us to compute LB sets using Algorithm 3. In column generation (see, e.g. Desrosiers and Lübbecke 2005), in each iteration a subset of promising columns is generated and appended to the restricted set of columns and the single objective linear program is resolved on . Column generation continues as long as new promising columns exist. Otherwise, an optimal solution has been found (i.e. in that case, an optimal solution of the single objective linear program on is also an optimal solution of the single objective linear program on ). Promising columns are identified using dual information. Let denote the dual variable associated with constraint (4) for a given and the dual variable associated with constraint (5), the pricing subproblem we have to solve corresponds to:
(9) 
subject to constraints (27)–(34) given in Appendix D, omitting subscript . It is an elementary shortest path problem with resource constraints that can be solved by means of a labeling algorithm (cf. Feillet et al. 2004). In our labeling algorithm, a label carries the following information: the node the label is associated with, the time consumption until that node, the reduced cost so far, which nodes have been visited along the path leading to the node, and a pointer to the parent label. In order to compute the reduced cost of the path associated with a given label, we use a reduced cost matrix that is generated before the labeling algorithm is called. The reduced cost of arc is , where if , is replaced by and gives the costs for traversing arc and the score for visiting location . Since the aim of this paper is not to investigate the most efficient pricing algorithm, we refrain from adopting all enhancements proposed in the literature (e.g. Righini and Salani 2009).
Note that during the execution of BIOBAB we keep all previously generated columns in the column pool. We only temporarily deactivate those columns that are incompatible with current branching decisions. When branching on objective space, it can happen that the current pool of columns does not allow to produce a feasible solution. This can be due to two reasons: either (i) because there exists no feasible solution to the current problem or (ii) because there exists a feasible solution but it cannot be reached with the currently available columns. In order to fix this issue and guarantee feasibility, we use a dummy column which allows to satisfy every constraint from the current problem, including the branching decisions described in Section 6.5.2. These include branching on control points. Therefore, some control points may be mandatory and others may be forbidden. Thus, the dummy column covers all mandatory control points, does not cover any forbidden control point, uses all the available vehicles, has a cost inferior to the maximum allowed cost and a profit superior to the minimum allowed profit. Using this column is penalized in the objective function, so that the column generation procedure converges to feasible solutions that do not use the dummy column. Assuming the variable associated to this dummy column is , the objective function described in Equation (8) is modified as follows:
(10) 
where is an arbitrarily large number. The dummy column is only activated when no feasible solution can be found, and it is systematically deactivated after column generation converges. At this point, if has a strictly positive value then there does not exist a feasible solution that satisfies all the branching decisions.
6.5.2 Branching scheme
Objective space branching is performed as usual (see Section 4.4). Since objective space branching involves setting bounds on both objectives, we include these constraints right from the start (initially setting the respective bounds to infinity) and later update them according to the branching decisions:
(11)  
(12) 
where is the local nadir point. In order to properly incorporate dual information from these two constraints into the subproblem, we modify the reduced cost matrix. Let be the dual variable value associated with constraint (11) and the value associated with constraint (12), then the reduced cost of a given arc is given by . Regarding decisionspace branching, we either branch on control points or on arcs (Line 4 in Algorithm 2). First, we check if a control point is visited a fractional number of times. Since each lower bound represents a set of solutions, the number of times a control point is visited may take different values within the same LB set. In order to select a control point to branch on, we consider the supported (fractional) solutions defining the current LB set. We then average the number of visits for each control point over this set of solutions. The control point with the average number of visits closest to is then selected for branching. In the case where each control point is visited 0 or 1 time on average, we check for arcs that are traversed a fractional number of times on average, following the same procedure. Branching on control points can be achieved by modifying constraint (4) in the master problem. In order to force the visit of control point , this constraint becomes . In order to forbid the visit of , this constraint becomes . As is usual in branchandprice for routing, branching on arcs involves adding constraints to the subproblem.
6.5.3 Applying the constraint method
In order to asses the quality of our BIOBAB on the BITOPTW, we compare it to the constraint method, using the faster direction. In this case, the function (Line 4 of Algorithm 4 in the Appendix) uses the same tree search algorithm as BIOBAB, i.e. the same rules for branching are applied, with the exception of objective space branching, since a singleobjective problem is solved. The bounding procedure is similar to the one used in BIOBAB, i.e. the same column generation code is used. However, instead of solving a succession of LPs to compute a LB set, a single solution is produced in each call to the procedure. Other than that, the code is exactly the same for both methods. This means that in the constraint framework, previously generated columns are reused. Additionally, a singleobjective solution to the linear relaxation yields a LB segment ; therefore, the existing UB set can be used to fathom subtrees using the filtering procedure described in Section 4.3. This works in combination with the fact that all produced integer solutions are stored, thus enhancing the constraint implementation with improvements developed for BIOBAB.
6.5.4 BITOPTW Benchmark instances
We use the TOPTW instances by Righini and Salani (2009) and we reduce their size by considering the first 15, 20, 25, 30 and 35 customers only. Using instances c101_100, r101_100, rc101_100 and pr01 and the number of customers mentioned above, we generate test instances for 1, 2, 3 and 4 vehicles. In total there are 80 BITOPTW instances.
6.5.5 Computational results
We now compare BIOBAB with the constraint method. Column generation can be timeintensive, so in this case all procedures are implemented in C++. The same column generation code is used by both BIOBAB and constraint. We run both algorithms on the same computer as mentioned above, and produce the performance profile depicted in Figure 16.
As we can see, BIOBAB offers a better performance overall. Although the constraint method is within the same order of magnitude for almost 90% of the instances, which is good, it can get up to 20 times slower in some rare cases. More strikingly, there are instances solved by BIOBAB that cannot be solved within the allotted CPU budget by the constraint method.
7 Conclusion
We have introduced a new algorithm for biobjective integer optimization. The algorithm is a generalization of singleobjective branchandbound to the biobjective context. Our algorithm is generic and uses a branching rule exploiting the biobjective nature of the underlying optimization problem. We have also developed several enhancements based on the integrality of the objective values of efficient solutions. Experiments show that these enhancements significantly improve the performance of BIOBAB, and make it competitive with stateoftheart generic methods for biobjective optimization for integer programs. In order to demonstrate the versatility of our framework, we have also applied it to a routing problem where the lower bound linear program is solved by means of column generation, thus providing, to the best of our knowledge, the first biobjective branchandprice algorithm. Standard criterion space search methods rely on the efficiency of MIP solvers and show their limitations in such cases: the proposed branchandbound algorithm outperforms the constraint method, which is a standard benchmark for multiobjective optimization.
One research perspective lies in generalizing the concept of lower bound set to more than two objectives, thus allowing exact approaches relying on this notion for more than two objectives. However this will raise the issue of the cardinality of the efficient set, which typically increases exponentially with the number of objectives. Therefore another research perspective lies in deriving approximation methods for multiobjective optimization based on exact procedures, with the goal of sacrificing neither the quality of solutions nor the readability of the produced set of solutions.
Acknowledgments
We wish two thank four anonymous referees for their valuable comments. Furthermore, financial support from the Austrian Science Fund (FWF): P23589N13 is gratefully acknowledged.
Appendix A Algorithms for criterion space search methods
Algorithms 4 and 5 outline the constraint method and the balanced box method, respectively. The set denotes the set of Pareto optimal solutions. Function can either be a call to Gurobi or CPLEX or to a tailormade branchandbound algorithm. The value of is supposed to be valid for the problem at hand; for instance with pure integer problems with integer coefficients for the objective function, is such a valid value.
Appendix B UBOFLP: Model formulation
In the UBOFLP, we are given a set of potential facilities and a set of locations. Furthermore, we denote by the set of locations that can be covered by facility because they are within a certain radius. Each facility has opening costs and each location has a weight or demand . The considered objectives simultaneously minimize the opening costs of facilities and maximize the total covered demand.
Using binary decision variables equal to if facility is opened and otherwise, and equal to if location is assigned to facility , we formally define the UBOFLP as follows:
(13)  
(14) 
subject to:
(15)  
(16)  
(17)  
(18) 
Appendix C SSUFLP: Model formulation
In the SSUFLP, we are given a set of potential facilities and a set of locations that have to be assigned to one of the facilities each. Each facility has opening costs and it costs to assign location to facility . The considered objectives simultaneously minimize the total opening costs of facilities and the total assignment costs.
Using binary decision variables equal to if facility is opened and otherwise, and equal to if location is assigned to facility , we formally define the SSUFLP as follows:
(19)  
(20) 
subject to:
(21)  
(22)  
(23)  
(24) 
Appendix D BITOPTW: Model formulation
The BITOPTW is defined on a directed graph , where is the set of arcs and the set of vertices, representing the starting location (vertex ), the ending location (vertex ) and control points. Each control point is associated with a score , a service time and a time window . Each route , with , has to start at location and end at location and each arc is associated with travel cost and travel time . The aim is to maximize the total collected score and to simultaneously minimize the total travel cost. Using binary decision variables equal to if location is visited and otherwise, and equal to if arc is traversed by route and otherwise, and continuous variables , denoting the beginning of service at by route , we formally define the BITOPTW as follows:
(25)  
(26) 
subject to:
(27)  
(28)  
(29)  
(30)  
(31)  
(32)  
(33)  
(34) 
Objective function (25) minimizes the total routing costs while objective function (26) maximizes the total collected profit. Constraints (27) and (28) make sure that each route starts at the defined starting point and ends at the correct ending point. Constraints (29) link the binary decision variables and (30) ensure connectivity for visited nodes. Constraints (31) set the time variables and (32) make sure that time windows are respected. We note that constraints (31) are not linear but they can easily be linearized using big terms.
References
 Adelgren et al. (2014) Adelgren, N., Belotti, P., and Gupte, A. (2014). Efficient storage of pareto points in biobjective mixed integer programming. arXiv:1411.6538 [cs.DS].
 Adelgren and Gupte (2016) Adelgren, N. and Gupte, A. (2016). A branchandbound method for multiobjective mixed integer linear programs and its implementation for biobjective problems. Working paper, Clemson University.
 Aneja and Nair (1979) Aneja, Y. P. and Nair, K. P. K. (1979). Bicriteria transportation problem. Management Science, 25:73–78.
 Baldacci et al. (2012) Baldacci, R., Mingozzi, A., and Roberti, R. (2012). Recent exact algorithms for solving the vehicle routing problem under capacity and time window constraints. European Journal of Operational Research, 218(1):1 – 6.
 Belotti et al. (2013) Belotti, P., Soylu, B., and Wiecek, M. M. (2013). A branchandbound algorithm for biobjective mixedinteger programs. available online at http://www.optimizationonline.org/DB_HTML/2013/01/3719.html.
 Belotti et al. (2016) Belotti, P., Soylu, B., and Wiecek, M. M. (2016). Fathoming rules for biobjective mixed integer linear programs: Review and extensions. Discrete Optimization, 22:341–363.
 Boland et al. (2015a) Boland, N., Charkhgard, H., and Savelsbergh, M. (2015a). A criterion space search algorithm for biobjective integer programming: The balanced box method. INFORMS Journal on Computing, 27:735–754.
 Boland et al. (2015b) Boland, N., Charkhgard, H., and Savelsbergh, M. (2015b). A criterion space search algorithm for biobjective mixed integer programming: The triangle splitting method. INFORMS Journal on Computing, 27(4):597–618.
 Braekers et al. (2016) Braekers, K., Hartl, R. F., Parragh, S. N., and Tricoire, F. (2016). A biobjective home care scheduling problem: Analyzing the tradeoff between costs and client inconvenience. European Journal of Operational Research, 248:428–443.
 Desrosiers and Lübbecke (2005) Desrosiers, J. and Lübbecke, M. E. (2005). A primer in column generation. In Desaulniers, G., Desrosiers, J., and Solomon, M. M., editors, Column Generation, pages 1–32. Springer US.
 Dolan and Moré (2002) Dolan, E. D. and Moré, J. J. (2002). Benchmarking optimization software with performance profiles. Mathematical Programming, 91(2):201–213.
 Ehrgott (2005) Ehrgott, M. (2005). Multicriteria optimization, volume 2. Springer.

Ehrgott and Gandibleux (2006)
Ehrgott, M. and Gandibleux, X. (2006).
Bound sets for biobjective combinatorial optimization problems.
Computers & Operations Research, 34:2674–2694.  Feillet et al. (2004) Feillet, D., Dejax, P., Gendreau, M., and Gueguen, C. (2004). An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems. Networks, 44(3):216–229.
 Gadegaard et al. (2016a) Gadegaard, S., Ehrgott, M., and Nielsen, L. (2016a). Biobjective branchandcut algorithms: Applications to the single source capacitated facility location problem. Working paper, Aarhus Universitet.
 Gadegaard et al. (2016b) Gadegaard, S., Nielsen, L., and Ehrgott, M. (2016b). An instance generator for the capacitated facility location problem. github, source code (v1.0.0).
 Haimes et al. (1971) Haimes, Y., Lasdon, L., and Wismer, D. (1971). On a bicriterion formulation of the problems of integrated system identification and system optimization. IEEE Transactions on systems, man and cybernetics, 1:296–297.
 Jozefowiez et al. (2012) Jozefowiez, N., Laporte, G., and Semet, F. (2012). A generic branchandcut algorithm for multiobjective optimization problems: application to the multilabel traveling salesman problem. INFORMS Journal on Computing, 24:554–564.
 Leitner et al. (2014) Leitner, M., Ljubić, I., and Sinnl, M. (2014). A computational study of exact approaches for the biobjective prizecollecting steiner tree problem. INFORMS Journal on Computing, 27(1):118–134.
 Masin and Bukchin (2008) Masin, M. and Bukchin, Y. (2008). Diversity maximization approach for multiobjective optimization. Operations Research, 56(2):411–424.
 Mavrotas and Diakoulaki (1998) Mavrotas, G. and Diakoulaki, D. (1998). A branch and bound algorithm for mixed zeroone multiple objective linear programming. European Journal of Operational Research, 107:530–541.
 Raith et al. (2012) Raith, A., Moradi, S., Ehrgott, M., and Stiglmayr,