{'fraction_non_alphanumeric': 0.041099432920183976, 'fraction_numerical': 0.00999316548755934, 'mean_word_length': 4.294152161157833, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An outer approximation algorithm for multi-objective mixed-integer linear and non-linear programming

In this paper, we present the first outer approximation algorithm for multi-objective mixed-integer linear programming problems with any number of objectives. The algorithm also works for certain classes of non-linear programming problems. It produces the non-dominated extreme points as well as the facets of the convex hull of these points. The algorithm relies on an oracle which solves single-objective weighted-sum problems and we show that the required number of oracle calls is polynomial in the number of facets of the convex hull of the non-dominated extreme points in the case of multiobjective mixed-integer programming (MOMILP). Thus, for MOMILP problems for which the weighted-sum problem is solvable in polynomial time, the facets can be computed with incremental-polynomial delay. From a practical perspective, the algorithm starts from a valid lower bound set for the non-dominated extreme points and iteratively improves it. Therefore it can be used in multi-objective branch-and-bound algorithms and still provide a valid bound set at any stage, even if interrupted before converging. Moreover, the oracle produces Pareto optimal solutions, which makes the algorithm also attractive from the primal side in a multi-objective branch-and-bound context. Finally, the oracle can also be called with any relaxation of the primal problem, and the obtained points and facets still provide a valid lower bound set. A computational study on a set of benchmark instances from the literature and new non-linear multi-objective instances is provided.


Introduction and motivation
Many practical problems involve several, often conflicting objectives, such as profitability or cost versus environmental concerns (Demir et al, 2014;Ramos et al, 2014;Eskandarpour et al, 2021) or customer satisfaction (Braekers et al, 2016). This means that there is, in the general case, no single optimal solution which optimizes all objectives simultaneously but a set of trade-off solutions which are better than all other solutions but incomparable among each other. A solution belonging to this set is called Pareto optimal and has the property that no objective function value can be improved without deteriorating another. The images of these solutions in the space of objective function values (the objective space or criterion space) are called non-dominated points and they together form the non-dominated frontier or Pareto frontier. In multi-objective optimization, we usually aim at identifying at least one Pareto optimal solution for each non-dominated point.
Our research aims to further advance general purpose exact methods for bi-and multi-objective (mixed) integer linear programming relying on the branch-and-bound (B&B) idea. All recent successful implementations rely on lower bound sets (Ehrgott and Gandibleux, 2007), see, e.g., Gadegaard et al (2019); Parragh and Tricoire (2019); Forget et al (2022); Adelgren and Gupte (2021) but they are currently either restricted to two objectives or they only address the pure integer and not the mixed integer case. Lower bound sets are sets in the objective space. Depending on the definition, they either contain the non-dominated frontier or a set of points which, together, dominate the non-dominated frontier. They are a natural multiobjective extension of the lower bounds obtained by, e.g., relaxations in single-objective optimization. In a bi-objective context, they are usually computed "from the inside" with so-called inner approximation (IA) schemes. Only Forget et al (2020Forget et al ( , 2022 use an outer approximation (OA) scheme to compute bound sets "from the outside". The bound sets in these two works are obtained from the linear programming (LP) relaxation of the MOMILP using an implementation of Benson's OA algorithm for multi-objective linear programming (MOLP).
The OA algorithm we are proposing can directly obtain such a lower bound set for a MOMILP and also some classes of multi-objective (mixed) integer non-linear programming problems without the need of a convex relaxation of the considered problem. The algorithm produces the non-dominated extreme points and the facets of the convex hull of these points for any number of objectives.
The study of efficient algorithms for the computation of extreme points and facets of multi-objective optimization problems is also of independent theoretical interest. It is well known that even for structurally simple MOLP problems, the number of extreme points can be exponential in the size of the input (Ruhe, 1988). The same observation is also true for many well-known multi-objective combinatorial optimization (MOCO) problems, e.g., the multi-objective versions of the assignment and shortest path problems. Consequently, a polynomial time algorithm for the problem classes our algorithm can tackle cannot exist. In recent years, a new complexity measure emerged that is capable of separating some of the complexities of these problem classes. Bökler (2018) shows that in an output-sensitive complexity measure, MOLP is an easy problem, while the more general MOMILP remains hard. The idea is to regard the running time not only as a function of the input, but also of the output size. This is a very useful measure for multi-objective optimization if the whole set of extreme points, facets, or the whole non-dominated set is sought. A central theorem in this field of research is that the extreme points of a MOCO problem can be computed efficiently, if the weighted-sum scalarization can be computed efficiently, i.e., in polynomial time (Bökler and Mutzel, 2015). The proof of this theorem roots in the study of MOLP, in which the efficient computation of non-dominated extreme points is possible. While there are efficient algorithms for the computation of non-dominated facets of an MOLP (Bökler, 2018), no such theorem could yet be proven for the facet computation of more general problem classes. Our work closes this research gap.

Problem definition
We are interested in general multi-objective optimization problems (MOP) of the following form and require further assumptions below: where f : X → R p , p ∈ N, and X is a set. For this problem we define Q + := cl convf (X ) + R p ≥0 , the Edgeworth-Pareto hull. We observe that Q + is a closed convex polyhedron. In our problem setting, we are interested in the facets of Q + . We make the following assumptions about the problem: 1. The value inf w T f (x) : x ∈ X is computable for every w ∈ Q p ≥0 and has finite value.
From its definition and (1) it follows that Q + is a rational polyhedron, rec Q + = R p ≥0 , and there is y ∈ R p with Q + ⊆ y + R p ≥0 , i.e., an ideal point exists. The assumptions do not generalize, i.e., there are MOPs that do not conform to these. Examples include multi-objective linear programming instances without an ideal point (conflict with (1)), but also continuous non-linear programming instances where Q + exhibits curvature (conflict with (2)).
On the other hand, the assumptions include MOMILP iff an ideal point is present, and also non-linear problems, e.g., multi-objective combinatorial problems with quadratic constraints and objective functions. Furthermore, non-convex problems can be solved in which curvature is present in the non-dominated set but not in Q + .

Contribution
In this work, we propose the first OA algorithm which computes the extreme points and facets of Q + for the problem defined in Section 1.1. Our algorithm is motivated by the OA approach of Benson (1998) for MOLP. It uses an oracle which consists of solving single-objective weighted-sum problems. We show that our algorithm needs a number of oracle calls which is polynomial in the number of facets of Q + for a large class of problems.
As a consequence for MOMILP, we are able to propose a similar theorem to the one by Bökler and Mutzel (2015): If the weighted-sum problem of a given MOMILP is solvable in polynomial time, the facets can be computed with incremental-polynomial delay. These results extend and complement the results of Bökler and Mutzel (2015), who showed that with an IA algorithm, the extreme points can be found with incremental-polynomial delay in this case.
From a practical perspective, next to providing a lower bound set at any point of its execution, each oracle call may also produce a new Pareto optimal solution. This makes its use within multi-objective B&B attractive also from the primal side. Finally, the algorithm also provides a valid lower bound set when the oracle is called with any relaxation of the original problem, such as the LP-relaxation (potentially augmented with valid inequalities) in case of MOMILP. This opens perspectives for multi-objective branchand-cut algorithms.
We observe that the problem we study also corresponds to the problem addressed in parametric integer linear programming with parametrization in the objective function (Geoffrion and Nauss, 1977), which is an interesting problem on its own. Moreover, a related problem is to compute the extreme points and/or facets of the convex hull of a linear projection of a mixed-integer problem onto a small subspace. Our algorithm is also capable of computing this projection. However, our focus in this paper will be on the optimization side.

Outline
In Section 2 we provide notation, definitions and other preliminaries needed in the remainder of the paper. Section 3 provides an overview of existing inner and outer approximation schemes. Section 4 provides a general outline of OA algorithms for multi-objective optimization. To work for a concrete problem class such as MOLP or MOMILP, a point separation oracle specific to the problem class is needed. In Section 5 we propose two generic separation oracles which only depend on the availability of a black-box algorithm for the single-objective problem of the considered problem class. The theoretical runtime of the algorithm when using the presented oracles is discussed in Section 5.4. Section 6 contains a computational study on instances from literature and also new non-linear instances, and also discusses implementation details. We give empirical results on the numerical accuracy of our approach and provide a comparison with PolySCIP, which is an IA solver for MOMILP. Finally, Section 7 concludes the paper.
2 Notation, definitions and preliminaries 2.1 General definitions For x, y ∈ R n , the relations x = y, x ≤ y, and x < y apply component-wise. We define R n ≥0 := {x ∈ R n : x ≥ 0}, the non-negative orthant. Let A ∈ Q m×n , C ∈ Q p×n , b ∈ Q m and n = n 1 + n 2 . If X = {x ∈ Z n1 × R n2 : Ax ≥ b} and f (x) = Cx, then the MOP is called MOMILP. If n 1 = 0, the problem is called a MOLP and for n 2 = 0, it is called a multi-objective integer linear programming problem (MOILP).
We call the set Q := f (X ) the feasible set in decision space. A point y ∈ Q is non-dominated, iff there is noŷ ∈ Q withŷ ≤ y and y =ŷ. A feasible solution x ∈ X is called Pareto-optimal, iff f (x) is non-dominated. A feasible solution x is a supported Pareto-optimal solution, iff there is a w ∈ R p with w > 0, such that x is an optimal solution to the weighted-sum problem or weighted-sum scalarization: For every Pareto-optimal x, we call f (x) a supported non-dominated point. Note that compared to MOLP, in MOMILP there can be non-supported non-dominated points. The existence of such points justifies the need for multi-objective B&B algorithms, with which they can be found. This also justifies the fact that in MOLP the Edgeworth-Pareto hull usually has a different name: It is called the upper image of the MOLP instance.

Polyhedra and faces
Let P = {Ax ≥ b} for A ∈ R m×n , b ∈ R m , n, m ∈ N be a polyhedron. Given a ∈ R n and α ∈ R, the set H a α := {a T x = α} is called a hyperplane. The hyperplane H a α is supporting P in F , iff F := P ∩ H a α = ∅ and for all x ∈ P , we have a T x ≥ α. For every not necessarily supporting hyperplane H a α , the set F := P ∩ H a α is called a face of P . A face with affine dimension 0 is called an extreme point, one with affine dimension dim P − 1 is called a facet. Since for fully dimensional polyhedra, the hyperplane supporting a polyhedron in a facet is uniquely determined up to scaling, we often identify a facet with its defining (in)equality. A vector r ∈ R n is called a ray of P , iff for every x ∈ P and γ > 0 the point x+γr ∈ P . The set of all rays of P is called its recession cone or rec P . If rec P = ∅, then we call P a polytope. For the Edgeworth-Pareto hulls we observe in our model, we always have rec Q + = R p ≥0 .

Supported and extreme points
Supported non-dominated points can also be characterized by the Edgeworth-Pareto hull: A point y ∈ Q is a supported non-dominated point iff it lies on the boundary of Q + . We call the extreme points of Q + the non-dominated extreme points of the problem as well as the facets of Q + the non-dominated facets of the problem. In our problem setting, we are interested in the facets of Q + as they provide us with a lower bound set; as well as in Pareto-optimal solutions x such that f (x) is a non-dominated extreme point, as they provide us with primal solutions. The polyhedron is illustrated in Figure 1 and the non-dominated extreme points are the red points with the black circles around them. For more background on multi-objective optimization in general, we refer to, e.g., Chinchuluun and Pardalos (2007); Ehrgott (2005).

Running times and output-sensitive complexity
Since the number of facets and extreme points of the Edgeworth-Pareto hull (also called parametric complexity) of many problems can grow exponentially, the traditional aim of a polynomial running time for algorithms computing this entities is not useful. Instead, we usually regard an algorithm as efficient, if its running time can be bounded by a polynomial in the input and the output size. We call such an algorithm an output-polynomial time algorithm. y 1 y 2 Figure 1: Exemplary illustration of the polytope Q + = cl conv(Q) + R 2 ≥ A more restrictive notion is an incremental-polynomial delay algorithm: The kth-delay of an enumeration algorithm is defined as the time between producing the kth and (k + 1)th output, including the time to produce the first output and the time after the last output until termination. In the incremental-polynomial delay model, we require the kth delay to be bounded by a polynomial in the input size and k. In contrast to output-polynomial time algorithms, incremental-polynomial delay algorithms are not allowed to print all outputs at the end, but at certain time points there have to be updates. However, the model still allows the delays between the outputs to grow depending on the size of the output so far. This is remedied by the most restrictive model we discuss: A polynomial delay algorithm is an enumeration algorithm where every delay is bounded by a polynomial in the input size.
Naturally, a polynomial delay algorithm is an incremental-polynomial delay algorithm and an incrementalpolynomial delay algorithm is also an output-polynomial algorithm but the reverse implications do not hold in general. See also Bökler et al (2017) for a more in-depth exposition of the subject.

State-of-the-art
In this section, we provide an overview on existing IA and OA algorithms for multi-objective optimization problems. A focus is on the state-of-the-art concerning finding the Edgeworth-Pareto hull of an MOMILP as the best achievable linear lower bound set for an MOMILP is the boundary of Q + . At the end of the section, we also discuss the existing research on bound computation for multi-objective (mixed) integer non-linear programming.
Since Q + is a polyhedron, there are two general ways to represent it: By its generators, i.e., extreme points and rays, or by inequalities. In the first representation, it suffices to compute the extreme points, since as we assume the feasible set to be bounded, the extreme rays of Q + are always the unit vectors. In a multi-objective or parametric optimization setting, the extreme point representation is more compelling, since it provides us with value vectors of efficient solutions. In our setting however, we are more interested in computing inequality representations. In bi-objective problems, there is not much of a difference between these representations, but with a higher number of objectives, the computational cost of switching between them grows exponentially.
To the best of our knowledge, the first works that were concerned with computing extreme points of Q + of MOMILP were independently conducted by Aneja and Nair (1979); Cohon (1978), and Dial (1979). Albeit, this algorithm, that is usually referenced as the dichotomic approach today, does only work in the bi-objective case. As the approach is well-known, we forgo describing it in detail. Although the approach has been formulated for transportation problems, it is indeed a general purpose tool. Much later, Bökler et al (2017) prove that the dichotomic approach can be implemented with polynomial delay. Moreover, in the special case of MOLP, this can even be achieved with polynomial space in the input size using a slightly different algorithm also reported therein.
In the seminal work of Benson (1998), the first outer approximation algorithm for solving MOLP problems in objective space is given. Despite its name, it is an exact method. It forms the basis of the open source vector linear programming solver bensolve (Löhne and Weißing, 2017) and motivates our work. The key idea of Benson's algorithm is that Q + can be approximated from the "outside", by adding new facets in each step, until the approximation S is equal to Q + . Initially, the approximation S = y + R p ≥0 , where y is the ideal point, and an interior pointr ∈ intQ + are required. Then, in each step a point that lies on the intersection of the boundary of Q + and the line connecting one of the current vertices y * of S andr is determined. It is used to obtain a supporting hyperplane of Q + , by solving an LP, or to prove that y * is a vertex of Q + . Bökler (2018) proves that this algorithm is capable of computing the extreme points of an MOLP in output-polynomial time and the facets even with incremental-polynomial delay.
Extending the work by Benson, Ehrgott et al (2012) develop a dual algorithm based on the geometric duality theory of MOLP by Heyde and Löhne (2008). Instead of operating in the decision or objective space, it operates in a dual space of the objective space. Thus the dual algorithm essentially is a version of Benson's outer approximation algorithm with a very crucial caveat: Bökler and Mutzel (2015) show that instead of solving an involved LP to find separating hyperplanes, in this representation, it suffices to solve weighted sum LPs. In spite of its dual description, the dual algorithm can be described as operating in the objective space in the following way: The initial dual outer approximation is essentially a non-dominated extreme point of the problem with an attached non-negative orthant in the objective space, R p ≥0 . Due to duality theory, choosing an extreme point of the current dual outer approximation is the same as choosing a facet of the current inner approximation in the objective space. Finding a new supporting hyperplane in dual space is again equivalent to finding a new supported point in the objective space. In this view, the dual of Benson's algorithm is also the first instance of what we now call an inner approximation for MOLP. Regarding theoretical running time, Bökler (2018) shows that the dual algorithm computes the extreme points of an MOLP with incremental-polynomial delay and the facets in output-polynomial time. Interestingly, the outer approximation is thus more efficient in computing facets and the inner approximation is more efficient in computing extreme points.
Very recently, Csirmaz (2021) introduced a unifying perspective on inner and outer approximation schemes for MOLP, giving skeletal algorithms for both. Each scheme starts from an initial approximation S of Q + . Then, the inner (outer) approximation scheme relies on a plane (vertex) separating oracle determining if a vertex (facet) is already part of Q + . If yes, it is marked at final, otherwise, S is updated. Output to either scheme is Q + in double description format (i.e., vertices and facets of Q + ). Csirmaz (2021) shows that Benson's algorithm and its variants follow the outer approximation scheme whereas the dual variant of Ehrgott et al (2012) falls into the described inner approximation scheme. While Benson's algorithm and its variants address the multi-objective linear case, also some algorithms for obtaining the extreme supported points of MOMILP have been proposed in the past years. To the best of our knowledge, none of them approximates Q + from the outside and as such classifies as an OA algorithm. Przybylski et al (2010) propose a recursive scheme that resorts to solving bi-objective problems via dichotomic search derived from the weights associated with extreme points of a given facet of the weight-set polytope. The weight-set is the set of all eligible weights in a weighted-sum scalarization. Decomposing it into weight-set components, for which a given extreme supported image is optimal, results in a weight-set decomposition. Computing this decomposition motivates the development of the proposed scheme. Przybylski et al (2010) show that their method works for three objectives and in theory also for more. An enumerative procedure is due toÖzpeynirci and Köksalan (2010). It relies on what the authors call stages which are defined by p points. In each step, the normal weight vector to the hyperplane H defined by a given stage R is used in the weighted sum scalarization of the MOMILP to obtain a new point y. If y corresponds to one of the points of R, R is identified as a facet defining stage. Otherwise, the new point y is added to the set of non-dominated points and the list of to be explored stages is updated with y. A main drawback of the method is that it requires the generation of dummy points to initialize the search. Bökler and Mutzel (2015) propose IA algorithms for MOCO problems. Their work builds on the dual variant of Benson's algorithm, observing that the LPs that need to be solved in the MOLP case can be replaced by solving weighted sum problems of the original MOCO problem, showing that the geometric duality theory for MOLP by Heyde and Löhne (2008) can also be applied to non-convex problems if the interest lies in the Edgeworth-Pareto hull of the problem. To overcome the drawback of enumerating redundant hyperplanes, they introduce a lexicographic variant which guarantees that only extreme points of Q + are identified. First complexity results are also due to their work: Bökler and Mutzel (2015) show that for every fixed number of objectives and MOCO with a polynomially solvable weighted-sum scalarization, the non-dominated extreme points can be computed with incremental-polynomial delay. This result can be straightforwardly extended to MOMILP under the assumption that an ideal point exists: By the same arguments of the classic proof that (the decision version of) MOMILP ∈ NP, we see that the encoding lengths of optimal extreme point solutions of the integer hull remain polynomially bounded in the input size. Optimizing over the integer hull of an MOMILP thus is analogous to the combinatorial case.
Borndörfer et al (2016) develop a general purpose solver named PolySCIP for MOILP and MOLP with an arbitrary number of objectives as part of the constraint integer programming suite SCIP. To compute the extreme points of the MOILP and MOLP, they rely on a lifted weight space polyhedron. This polyhedron is essentially the dual polyhedron from geometric duality of MOLP applied to MOMILP, extending the work of Bökler and Mutzel (2015). PolySCIP employs a dual Benson outer approximation algorithm on this lifted weight space/dual space. Hence as discussed above, PolySCIP also falls in the category of inner approximations.
Even more recently, Halffmann et al (2020) have proposed an inner approximation algorithm for threeobjective mixed integer linear programs, with the purpose of determining the weight-set decomposition. The proposed algorithm relies on solving bi-objective subproblems via the dichotomic approach to identify these components. It runs in output polynomial time. Finally, building on the work ofÖzpeynirci and Köksalan (2010) and Przybylski et al (2010), Przybylski et al (2019) have proposed a new IA algorithm for MOMILP. In each step, it generates the convex hull of the supported extreme points which have been identified so far. Then, a yet unexplored facet is chosen and used to define the next weighted sum scalarization to be solved. The solution may either be a new non-dominated supported point or identify the facet to be (partially) part of Q + . The algorithm ends once all facets have been explored.
Multi-objective non-linear programming problems have received much less attention than linear problems. Klamroth et al (2003) develop inner and outer approximation algorithms producing piecewise linear approximations of the Pareto frontier for convex as well as nonconvex problems. An outer approximation algorithm for convex multi-objective programming problems with convex Pareto frontiers is due to Ehrgott et al (2011). It is an extension of Benson's OA algorithm and provides a set of weakly ǫ-non-dominated points. Building on this work, Löhne et al (2014) extend Benson's outer approximation algorithm and its dual variant to convex vector optimization. Niebling and Eichfelder (2019) use Benson cuts to improve the lower bound, which is initially based on the ideal point, in the context of a branch-and-bound scheme for smooth convex multi-objective optimization. More recently, an outer approximation based B&B algorithm has been proposed by De Santis et al (2020) for multi-objective mixed integer convex optimization problems. In the lower bound computation, building on the work of Ehrgott et al (2011) and Löhne et al (2014), a convex relaxation-based outer approximation is used. It relies on the computation of supporting hyperplanes. This computation is steered by the available local upper bounds within the B&B algorithm.

The outer approximation algorithm
In this section, we first give a general outline of OA algorithms for general multi-objective optimization problems. To this end, we generalize the terminology introduced by Csirmaz (2021). The key ingredient in implementing such algorithms for concrete problem classes like MOLP or MOMILP is the specification and implementation of a point separating oracle for the respective problem class.
Definition 1 (Point separating oracle (Csirmaz, 2021)). A point separating oracle for a polytope Q + ⊂ R p is a black box algorithm which takes as input a point y * ∈ R p and returns as tuple (status, H). The output is as follows: i) (inside, ∅), if y * ∈ Q + , or ii) (outside, H), where H is a supporting hyperplane H = {y ∈ R p : w T y = α} of Q + such that w T y * < α and w T y ≥ α for each y ∈ Q + (i.e., H separates y * from Q + ).
The OA algorithm is described in Algorithm 1. The approximations S = S 0 ⊃ S 1 ⊃ . . . used within the algorithm are considered to be stored in a double description format, i.e., as vertices and facets. The algorithm starts with an initial approximation S. This initial approximation consists of the ideal point y I ∈ Q p , whose ith component is the minimum value in the ith objective, together with the non-negative orthant, i.e., S := y I + R p ≥ . It proceeds in an iterative fashion by checking the vertices of the current approximation S i for containment in Q + using the separation oracle. This is done repeatedly until for some S i ′ the oracle answers with status inside for each vertex, which means S i ′ = Q + .
Algorithm 1: Generic OA algorithm for multi-objective optimization Data: initial approximation S specified by double description (vertices and facets) Result: Q + specified by double description Two steps of the OA algorithm are illustrated in Figure 2. It shows the image of all feasible solutions to an MOILP as green dots. Figure 2a depicts the initial approximation S 0 , consisting of the ideal point y I and R p ≥ , in green. Then, the point separating oracle is called and the hyperplane H is returned which separates y I from Q + , as shown in Figure 2b. The new approximation S 1 has two vertices depicted in blue vertices(S 1 ) = {y A , y B }. An oracle call for the vertex filled with green, y A , will return inside, since it belongs to Q + , whereas for the second vertex y B , the oracle will produce a new hyperplane, separating y B from Q + .
The convergence-behaviour of an OA algorithm depends on the separation oracle. In the original work by Benson (1998) for MOLP, the supporting hyperplanes to Q + were only face supporting. The first facet supporting hyperplane producing oracles were given in Bökler (2018); Csirmaz (2021) for MOLP. Hence, we are interested in facet supporting hyperplanes and we show in the following sections that our proposed general separation oracle produces facets of Q + .
Remark 2. It is easy to see, that at any point of the algorithm, we have Q + ⊆ S i . Moreover, the algorithm still gives an OA of Q + when the separation oracle does not separate with respect to Q + but to any superset S(Q + ) ⊃ Q + . This means that Q + can be replaced with any relaxation of the polyhedron Q + in the separation oracle when the goal is to obtain lower bound sets for multi-objective B&B algorithms.
Remark 3. Let us briefly discuss the key running-time insights of the general OA algorithm. Let us therefore assume, that a point separating oracle is available, produces facets only, and its running time is a constant as a simplification. In this case it is easy to see that for every extreme point found, i.e., the situation of the first if-branch, we ask the point separating oracle once. And for every facet found, i.e., the situation of the else-branch, we also ask the oracle once. No facet and no extreme point is found twice. Thus, we see that the overall running time is polynomial in the number of facets, extreme points, and the running time of the point separating oracle. Moreover, every time we find a facet it is clear that we can output it directly. This is not true for extreme points though, as the extreme points found in the if-branch are only intermediate extreme points and may be deleted in later iterations. The real extreme points of the Edgeworth-Pareto hull are only know at the time the algorithm terminates. This is the main reason why the algorithm is considered more efficient in computing facets than extreme points. Consequently, a finite number of extreme points and facets is needed unless we want to approximate the Edgeworth-Pareto hull. For an in-depth look into the running-times of the inner and outer approximation algorithms see Bökler (2018).

Geometric duality
The techniques we use in this paper are heavily influenced by the theory on MOLP. For MOLP, a duality theory (Heyde and Löhne, 2008) exists that is built around a polyhedron D, the lower image, that is geometrically dual to Q + . For w ∈ R p−1 , we define λ(w) := w 1 , . . . , w p−1 , 1 − p−1 i=1 w i ∈ R p . For an MOLP with A ∈ Q m×n , b ∈ Q m , and C ∈ Q p×n , the lower image is defined as D := (w 1 , . . . , w p−1 , α) ∈ R p : We can characterize D independently of A, b, and C and only with respect to the upper image polyhedron Q + in the following way: In this fashion, we give a more general definition of the lower image for any MOP with our additional requirements: Definition 4. Let Q + be the Edgeworth-Pareto hull of an instance of a MOP with the above requirements. The generalized lower image is defined as: Observe that every such finitely generated polyhedron Q + can be perceived as the upper image of some MOLP, e.g., the decision polyhedron is also Q + and the objective matrix is the identity. Thus, we can apply the geometric duality theory to general Edgeworth-Pareto hulls as defined above. In particular, there exists an inclusion reversing one-to-one map from the k-dimensional faces of Q + to the (p − k − 1)-dimensional faces of D. For our purpose, the most important corollary is the following: For every extreme point ( w 1 , . . . , w p−1 , α) of D, the hyperplane y ∈ R p : λ( w) T y = α is facet supporting to Q + .

A first general separation oracle
Let w i be the coefficient of objective i, i = 1, . . . , p in the desired hyperplane H and α be the right-hand-side of H.
Problem (Sep-y * ) encodes that (w, α) should induce a separating hyperplane by enforcing that all y ∈ Q should be on the positive side of the hyperplane using constraints (FEAS). Constraint (WSUM1) is a normalisation constraint for the coefficients of the obtained hyperplane. It will be used in the proofs later on. If the objective function value of (Sep-y * ) is negative for a given y * , we get that y * ∈ Q + and (w, α) gives the coefficients of the corresponding separating hyperplane. To deal with the infinite size of Q, we propose to solve (Sep-y * ) using a cutting-plane approach, where constraints (FEAS) get separated: Given a solution ( w, α), there exists a violated constraint (FEAS), iff inf w T y : y ∈ Q = inf w T f (x) : x ∈ X < α. Thus, the separation problem for (FEAS) consists in solving a weighted-sum problem.
Remark 5. In our implementation, the weighted-sum problems are weighted-sum scalarizations of the given original mixed-integer problems. We thus obtain integer feasible solutions. Further implementation details are given in Section 6.2.
We note that the separation oracle (Sep-y * ) shares some similarity with the LPs used in local cut separation for single-objective mixed-integer linear programming (Applegate et al, 2001;Chvátal et al, 2013).
We now prove that this is an implementation of a general point separating oracle and that new separating hyperplanes are facet supporting.
Proof. If (y * ) T w − α < 0, then w T (y * ) < α and w T y ≥ α for all y ∈ Q + by constraints (FEAS). In other words, H w, α := y ∈ R p : w T y = α separates y * from Q + and thus y * / ∈ Q. If y * is not in Q + , then by the hyperplane separation theorem, there is a hyperplane H w, α , separating y * from Q + , i.e., w T (y * ) < α and w T y > α for all y ∈ Q + . Moreover, we can assume w to be normalized, i.e., p i=1 w i = 1. As rec Q + = R p ≥0 , we can also assume all w i to be at least 0. Consequently, ( w, α) is a feasible solution to Sep-y * and its value (y * ) T w − α is at most 0 and an upper bound on the optimal value. Lemma 7. For every optimal extreme point solution ( w, α) of Sep-y * , we have H w, α := y ∈ R p : w T y = α is facet supporting for Q + .
Proof. By settingᾱ := (y * ) T w − α and using that rec Q + = R p ≥0 , we can rewrite the feasible set of Sep-y * as follows: (w 1 , . . . , w p ,ᾱ) : ∀y ∈ Q + : w T (y − y * ) +ᾱ ≥ 0, w i ≥ 0, finding the smallestᾱ in the above set is equivalent to finding the largestᾱ in (w 1 , . . . , w p−1 , −ᾱ) : ∀y ∈ Q + − y * : λ(w) T y ≥ᾱ, w i ≥ 0, which is the generalized lower image of Q + − y * . An optimal extreme point solution of Sep-y * is thus an extreme point of the corresponding generalized lower image. The claim thus follows from the geometric duality theorem.
Theorem 8. The above algorithm using (Sep-y * ) is a point separating oracle that provides us with facet supporting inequalities.

A target cut-like separation oracle
Another separation oracle can be defined as follows. This oracle follows the target-cut paradigm introduced in Buchheim et al (2008) for single-objective mixed-integer linear programming. In the description below, we assume that all points encountered in the separation process are ≥ 1 (component-wise). This is without loss of generality, as we can always add a large constant to all points. In Section 6.2 we describe how we do this in our implementation. In the presentation of this oracle we start from the dual side and then proceed to the primal side.
Constraints (D-T-FEAS) impose that y * must be equal to a point, or dominated by a point which can be obtained by a linear combination of the points in Q. Constraints (TSUM) ensure that this linear combination only has positive coefficients. The key insight is that y * only belongs to Q + iff the objective of (D-TSep-y * ) is one or larger, see Buchheim et al (2008) for more details.
(TSep-y * ) min (y * ) T w s.t. y T w ≥ 1 ∀y ∈ Q (TFEAS) Using the above ideas, we can implement a point separating oracle as follows. If an optimal solution w of (TSep-y * ) has an objective value smaller than 1, w T y ≥ 1 gives a separating hyperplane. Otherwise, we return inside. Similar to constraints (FEAS) of (Sep-y * ), we propose to solve (TSep-y * ) using a cutting-plane approach, where constraints (TFEAS) are separated by solving weighted-sum problems. The reason this algorithm is a point separating oracle is analogous to Lemma 6: We see that for any feasible solution w with objective value (y * ) T w < 1, the hyperplane {x ∈ R p : w T x = 1} again separates y * from Q, as long as Q ⊆ R p ≥ . Moreover, in case the minimum value is at least 1, then there is no separating hyperplane, hence y * ∈ Q + .
To prove that an optimal extreme point solution corresponds in fact to a facet supporting inequality, we observe that the feasible set of TSep-y * is a polar dual polyhedron to Q + .
Theorem 9. The above algorithm using (TSep-y * ) is a point separating oracle that provides us with facet supporting inequalities.

Theoretical running time
Let us first analyze the most general version of the OA algorithm. Using the methodology of Bökler (2018), we observe the following corollary: Corollary 10. The OA algorithm with the point separating oracles from Section 4 produces the next facet in time O k ⌊ p 2 ⌋ (T O + k log k) , where T O is the running time of the oracle and k ∈ N is the number of facets already computed.
Both our new point separating oracles show that a point separating oracle can be implemented by solving a sequence of weighted-sum scalarizations. The major advantage of solving weighted-sum scalarizations over other known methods to implement the oracle model is that no new constraints are added and thus the structural properties of the feasible sets remain untouched. It remains to show, however, that the number of weighted-sum scalarizations that need to be solved do not grow too fast. The problems Sep-y * and TSep-y * are linear programming problems and thus only weakly polynomial-time algorithms are known to solve them. Hence, we need to make some assumptions about the encoding length of their input, while we first avoid to restrict the running time of the weighted-sum oracle.
Lemma 11. An optimal extreme point of (Sep-y * ) and (TSep-y * ) can be computed with a number of weightedsum calls that grows polynomially in the encoding length of the extreme points of Q + if an upper bound ϕ ∈ N on the maximum of these encoding lengths is available.
Proof. Using the methodology of Grötschel et al (1993), we compute an optimal extreme point solution to the LPs Sep-y * or TSep-y * in polynomial time by using the ellipsoid method and employing a rounding mechanism to reach an optimal extreme point solution. The ellipsoid method has a running time polynomial in the number of variables (p+ 1) and the encoding length of extreme points of the input polyhedron, i.e., the feasible set of Sep-y * or TSep-y * . The encoding length of the extreme points of the feasible sets of Sep-y * and TSep-y * are polynomially bounded in the encoding length of their facets. The facets of Sep-y * and TSep-y * , in turn, correspond to the extreme points of Q + . The bound on the encoding length ϕ is therefore necessary for the ellipsoid method to bound the volume of the feasible set Sep-y * and TSep-y * and thus the number of its iterations.
The separation oracle in the invokation of the ellipsoid method is implemented by a solver for the weighted-sum problem. The weighted-sum solver is hence called only polynomially many times.
For a large problem class the extreme points of Q + can be bounded by a polynomial in the input size, namely MOMILPs: It is well-known that the encoding length of every extreme point of the integer hull is bounded by a polynomial in the input size. The encoding length of the image of each such feasible extreme point under the objective function is thus also polynomially bounded. In fact, the objective functions do not need to be linear; polynomial computability in extreme points of the integer hull suffices. In all these cases, ϕ in the above lemma can be estimated using these bounds.
Corollary 12. For MOMILP, an optimal extreme point of (Sep-y * ) and (TSep-y * ) can be computed with a polynomial number of weighted-sum calls.
We now proceed and add the restriction of a polynomial-time solvable weighted-sum problem. In this case, the whole point separation can be performed in polynomial time for any MOMILP and the following theorem follows from the above.
Theorem 13. If the weighted-sum problem for a given MOMILP is polynomial time solvable, then the facets of Q + can be computed with incremental-polynomial delay.
This theorem complements the work by Bökler and Mutzel (2015), where the same result was shown for linear combinatorial optimization problems and extreme point of Q + . It is easy to see, that the result can also be applied to MOMILPs.

Computational experiments
The proposed OA algorithms are implemented in Python using CPLEX 12.10 as LP/mixed integer (non-)linear programming-solver and the parma polyhedral library (ppl) (Bagnara et al, 2008) for vertex enumeration. In this section, we first describe the instances used in our computational experiments in Section 6.1, followed by a discussion of implementation details in Section 6.2 and results in Section 6.3.

Instances
We use instance sets from the literature and also new non-linear instances. Our newly created instances are available online under https://msinnl.github.io/pages/multi-objective.html. The sets are as follows.
• AP-MOO: These are multi-objective assignment instances from the multi-objective optimization library (moolibrary) 1 , with three objectives. The number of agents and tasks is identical and ranges from 5 to 40, in increments of 5. The objective function coefficients are random integers in the range [1,20].
There are 100 instances in the set. They were proposed in Kirlik and Sayın (2014); similar instances were also used in experiments for IA algorithms in (Özpeynirci and Köksalan, 2010;Przybylski et al, 2010Przybylski et al, , 2019.
• KP-MOO: These are multi-objective knapsack instances from the moolibrary with three to five objectives. Both profits and weights are random integers from the interval [1,1000]. The budget is calculated as half the total weight of all items, rounded up to the next integer. There are 160 instances in the set: 100 instances with three objectives (10 to 100 items), 40 with four (10 to 40 items), and 20 with five (10 to 20 items). Like the set AP-MOO, these instances are proposed in Kirlik and Sayın (2014) and similar instances are also used in experiments for IA algorithms in (Özpeynirci and Köksalan, 2010;Przybylski et al, 2010Przybylski et al, , 2019.
• QCP-MOO: These are three-objective quadratic covering instances newly created for this work. The problem is defined as follows: , 2, 3}, w ∈ Z n and C ∈ Z. We have created ten instances for each n ∈ {10, 20, 30, 40}. The entries of V j are random integers from the interval [1, 10] and the entries of w are random integer from the interval [1,100]. The value of C is set to ⌊w T x/4⌋. The structure of the problem allows an easy linearization of the non-linear objective function using a standard McCormicklinearization to transform the problem into an MOMILP. This allows us to analyze potential benefits of solving the problem directly as a non-linear multi-objective problem against solving the linearized version of it.

Implementation details
Instance transformation and offset calculation. We first transform the instances to minimization form by flipping the objective coefficients for all objectives which are maximization. We then compute y I by solving min y∈Q e T i y for unit-vector e i . If the value of y I is negative for an objective j, we add the constant −y I j + 1 to the respective objective to ensure that all points encountered in the separation process are component-wise ≥ 1 as needed for (TSep-y * ). Moreover, instead of using the value of one on the righthand-side of constraints (TWSUM), we fix the value to the sum of the components of the ideal point (after adding the constant −y I j + 1 to objectives j with negative value). This proved to be more numerically stable in preliminary computations.
Separation. For the initial approximation S, the ideal point y I can be obtained by solving min y∈Q e T i y for unit-vector e i . Moreover, we initialize the separation LPs by adding constraints (FEAS), resp., (TFEAS) induced by the solutions obtained for the p problems used for calculating the ideal point. When separating constraints (FEAS), resp., (TFEAS), we use tolerance ε for checking violation. Once a constraint (FEAS), resp., (TFEAS) is added to the separation LP, we leave it there for the remainder of the algorithm (i.e., for all subsequent separation oracle calls). In the separation LPs, the numerical emphasis parameter of CPLEX is turned on, and the feasiblity and optimality tolerances are set to 1e-9, i.e., the most accurate value possible. When checking the result of the separation oracle (i.e., the objective value of the separation LP against zero, resp., one), we also use tolerance ε.
We use ppl with integer numbers as input, as using fractional numbers leads to (more) numerical instabilities. In order to do so, we scale each (w, α) obtained from the separation oracle by 10 9 and take the integer part of the obtained number. When checking if a point y * ∈ vertices(S i ) is in Q + , we proceed as follows: When a point y * is within tolerance value ε to a point in y ′ ∈ insideV ertices for each coordinate, we consider y * = y ′ and do not call the separation oracle for it. For each S i , we check the points in the order in which they are made available by ppl. Moreover, when updating the outer approximations, we use a slightly different strategy compared to the outline in Algorithm 1: We do not immediately calculate S i ∩ H and move to S i+1 once we discovered a y * ∈ insideV ertices; instead, we check all y * ∈ vertices(S i ), collect the obtained separating hyperplanes, and then add them all at the same time to obtain the next outer approximation. This approach turned out to be faster and more numerically stable in preliminary computations. Finally, we also add every point found when separating (FEAS), resp., (TFEAS) to insideV ertices, as these points are obtained from the solution of a weighted-sum problem, and thus are supported non-dominated points.

Results
Numerical accuracy of our approach and the influence of ε. The focus in the implementation of our algorithms is to find a good balance between numerical accuracy and speed. We note that IA algorithms can suffer from numerical issues, for example in Przybylski et al (2019) different versions of the same IA algorithm are tested, and some versions cannot find the complete set of non-dominated extremes points for some instances. In our computational study, we use the open-source IA solver PolySCIP (Borndörfer et al, 2016;Maher et al, 2017) for comparison. As experiments show (see below), there are some instances were PolySCIP does not find all extreme points; moreover, we observed in PolySCIP's output that in some cases it also reports some weakly-dominated points as part of the set of non-dominated extremes points. In OA, numerical issues can lead to missing cut-off points, which do not belong to Q + . However, even if we miss some of these points, the obtained solution is still a valid outer approximation.
To investigate the influence of different values of ε on the performance and accuracy of our approach, we first consider the instance set AP-MOO. Since the assignment problem has a totally unimodular constraint matrix, it can be solved as a linear programming problem and thus we can obtain the extreme points of Q + by using the multi-objective LP solver bensolve. Table 1 shows the number of extreme points of Q + obtained by our OA with both oracles (columns |Q + |), and the number of extreme points when these points are rounded to the next integer (columns |Q + |, recall that AP-MOO has integer coefficients, thus all extreme points have integer values). We considered ε = 1e − 3 and ε = 1e − 5 for both oracles. We also report the number of extreme points obtained by bensolve and PolySCIP and the run time of the methods (column t[s]). The results are aggregated by instance size. To allow for a meaningful comparison, we only report results for instance sizes were all approaches terminate within the time limit for all instances of this size.
The results show that for instances with size up to 30, all methods give a consistent number of extreme points. For instances of size 35, the size of |Q + | obtained by the OA is consistent with bensolve considering both oracles when ε = 1e − 5. For instances of size 40, the number of extreme points obtained by the OA and by bensolve becomes inconsistent, but only by a small amount. This is not surprising as with larger instances size, the size of |Q + | also grows, and thus the potential for numerical inaccuracies increases as well. We observe that the number of non-rounded extreme points (|Q + |) also stays fairly consistent for smaller-sized instances, and for all sizes is about three times as large as the number of points after rounding. Regarding the run time, using ε = 1e − 3 is about twice as fast as using ε = 1e − 5.  Results for instance set KP-MOO. Table 2 shows the results for instance set KP-MOO obtained by our algorithms and by PolySCIP. Our algorithms are run with ε = 1e − 5. In these tables we give the number of rounded extreme points (columns |Q + |), the number of facets (columns #fac) and the runtime (columns t [s]). We see that the number of points is quite similar for all three methods, however, as the instances become larger for p = 3 and for p = 4, 5 in general, the numbers differ. This is not unexpected, as both larger instances and more objectives could have a negative effect on numerical accuracy. In particular, as for larger objectives, the number of facets is growing considerably. We observe that for the number of facets, the difference is larger compared to the number of points. This could be explained by the fact that |Q + | is the rounded set of points, which exploits that the instance coefficients are integer. Thus, numerical inaccuracies are corrected by the rounding, e.g., multiple slightly fractional points are likely to be rounded to the same integer point. The runtime-performance of both separation oracels is quite similar. Results for instance set QCP-MOO. Table 3 shows the results for instance set QCP-MOO obtained by our algorithms in both the variant where the problem is solved directly as a multi-objective mixed-integer nonlinear programming problem, and also the variant where the problem is linearized, i.e., where a MOMILP gets solved. This means that the black-box solver for the weighted-sum problems must either solve mixed integer non-linear programming problems or mixed integer linear programming problems. We note that CPLEX is capable of both. We report the same values as in Table 2 and additionally report the number of instances solved out of ten for each item-size (columns #sol). We see that for these instances, the oracle Sep-y * seems more efficient, as both in the non-linear as well as in the linearized case two more instances can be solved. By looking at the item-sizes where all instances can be solved (i.e., 10 and 20) the numerical accuracy does not seem to be affected by the type of black-box solver used. On the other hand, the results show directly solving the non-linear programming variant of the problem instead of the linearized is generally faster. This means it could pay off that our proposed algorithm is capable of directly solving multi-objective mixed-integer non-linear programming problems without the need for linearization.

Conclusion
We present the first outer-approximation algorithm to compute the extreme points and facets for multiobjective mixed-integer linear programming problems and also certain multi-objective mixed-integer non-linear programming problems. We show that the number of weighted-sum mixed-integer oracles needed to compute these extreme points and facets can be bounded polynomially by the size of the output and the input for MOMILP. And for MOMILP with polynomial-time computable weighted-sum scalarizations, the facets of the Edgeworth-Pareto hull can be computed with incremental polynomial delay. We provide a computational study on instances from literature and new non-linear instances. In this study, we investigate the numerical accuracy of our method and also include a comparison with PolySCIP, which is an inner approximation algorithm for multi-objective integer linear programming problems. We conjecture that existing and future multi-objective B&B methods can benefit from partial bound set computation which now is possible with our proposed method. Future work should investigate the efficacy of this approach in practical B&B algorithms.
Data availability The datasets generated during and/or analyzed during the current study are available in the web pages http://home.ku.edu.tr/~moolibrary/ and https://msinnl.github.io/pages/multi-objective.html.