{'fraction_non_alphanumeric': 0.0461954197951051, 'fraction_numerical': 0.012468145513006892, 'mean_word_length': 4.079047949965253, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Outer approximation algorithms for convex vector optimization problems

In this study, we present a general framework of outer approximation algorithms to solve convex vector optimization problems, in which the Pascoletti-Serafini (PS) scalarization is solved iteratively. This scalarization finds the minimum 'distance' from a reference point, which is usually taken as a vertex of the current outer approximation, to the upper image through a given direction. We propose efficient methods to select the parameters (the reference point and direction vector) of the PS scalarization and analyze the effects of these on the overall performance of the algorithm. Different from the existing vertex selection rules from the literature, the proposed methods do not require solving additional single-objective optimization problems. Using some test problems, we conduct an extensive computational study where three different measures are set as the stopping criteria: the approximation error, the runtime, and the cardinality of solution set. We observe that the proposed variants have satisfactory results especially in terms of runtime compared to the existing variants from the literature.


Introduction
Multiobjective optimization refers to optimizing multiple conflicting objectives simultaneously and it is a useful tool for many applications in various fields from finance to engineering since, by the nature of applications, there may be a trade-off among several objectives. For a multiobjective optimization problem (MOP), there is no single solution optimizing all objectives. Instead, feasible solutions that cannot be improved in one objective without deteriorating in at least another objective, namely efficient solutions, are of interest. In many applications, instead of finding the set of all efficient solutions, it is sufficient to obtain the images of efficient solutions, namely the nondominated points, in the objective space.
An optimization problem that requires minimizing or maximizing a vector-valued objective function with respect to a partial order induced by an ordering cone C is referred to as a vector optimization problem (VOP) and it is also widely used in different fields: for instance, see [11,22] for applications in financial mathematics. MOPs can be seen as special instances of VOPs where the ordering cone is the nonnegative orthant. The terminology in vector optimization is slightly different from the multiobjective optimization terminology. In particular, for a minimization problem where the ordering cone is C, a minimizer is a feasible solution that cannot be improved with respect to the order relation induced by C. The image of a minimizer is then a C-minimal point in the objective space.
One of the most common approaches to obtain a minimizer for a VOP is to solve a scalarization problem, which is a single objective optimization problem induced by the VOP. In general, a scalarization model is parametric and has the potential to generate a 'representative' set of minimizers when solved for different sets of parameters. Throughout, two scalarization methods will be used. The first one is the well-known weighted sum scalarization [8], which is performed by optimizing the weighted sum of the objectives over the original feasible region. In addition, we use the Pascoletti-Serafini (PS) scalarization [20], which aims to find the closest C-minimal point from a given reference point through a given direction parameter. Unlike weighted sum scalarization, it has the potential to find all minimizers of a given VOP even if the problem is not convex.
In this paper, we focus on convex vector optimization problems (CVOPs). In the literature, there are iterative algorithms utilizing scalarization models to generate an approximation of the set of all C-minimal points in the objective space. To the best of our knowledge, the first such algorithm is proposed in 1998, by Benson, [4] and it is designed to solve linear MOPs. It generates the set of all nondominated points in the objective space by iteratively obtaining improved polyhedral outer approximations of it. Later, this algorithm is generalized to solve linear VOPs; moreover, using geometric duality results, a geometric dual counterpart of this algorithm is also proposed, see [15]. In 2011, Ehrgott et al. [7] extended the linear VOP algorithm from [15] to solve CVOPs. Then, in 2014, Löhne et al. [16] proposed a similar algorithm and a geometric dual variant.
On the other hand, in 2003, Klamroth et al. [14] proposed approximation algorithms for convex and non-convex MOPs. The mechanism of the outer approximation algorithm for the convex case is similar to Benson's algorithm [4]. Different from [4], in [14] the selection for the reference point for the scalarization is not arbitrary. In [14], there is also an inner approximation algorithm for the convex problems and the convergence rates of both algorithms are provided for the biobjective case.
Recently, Dörfler et al. [6] proposed a variant of Benson's algorithm for CVOPs that includes a vertex selection procedure that also yields a direction parameter for the PS scalarization. To compute the parameters, the algorithm solves quadratic programming problems for the vertices of the current outer approximation.
In addition to the algorithms which solve PS scalarization (or equivalent models), recently Ararat et al. [2] proposed an outer approximation algorithm that solves norm-minimizing scalarizations. This scalarization does not require a direction parameter but only a reference point, which is again selected among the vertices of the current outer approximation.
The aforementioned algorithms differ in terms of the selection procedures for the parameters of the scalarization. The ones from [14,6] require solving additional models to select a vertex at each iteration. This feature enables them to provide the current approximation error at each iteration of the algorithm at the cost of solving a considerable number of models. The rest of the algorithms do not solve additional models and they provide the approximation error after termination. For the selection of the direction parameter, [14,7] use a fixed point from the upper image; [6] uses the point from the inner approximation that yields the minimum distance to the selected vertex; and [16] uses a fixed direction from the interior of the ordering cone through the algorithm. Table 1   The main contributions of this paper can be listed as follows: (1) We present a general framework of the outer approximation algorithms from the literature; (2) propose various additional variants, which select the two parameters of the PS scalarizations in structured and efficient ways; and (3) compare the performances through numerical tests. In particular, after proposing different direction selection rules and conducting a preliminary computational study, we propose three vertex selection rules. Different from the vertex selection rules from the literature, these rules are not based on solving additional single-objective optimization problems, hence they are computationally more efficient. The first one benefits from the clustering of the vertices. The second rule selects the vertices using the adjacency information among them. For the last one, we employ a procedure that creates local upper bounds to the nondominated points. This procedure is first introduced in the context of multiobjective combinatorial optimization, see for instance [5], and by its design, it works only for convex MOPs. Together with the proposed variants, we also implement some of the algorithms from the literature and we provide an extensive computational study to observe their behavior under different stopping conditions. In Sections 2 and 3, we provide the preliminaries and the solution concepts together with the scalarization models used in the paper, respectively. In Section 4, we present the general framework of the outer approximation algorithm; the construction of different direction and vertex selection rules; and similar algorithms from the literature. We provide the test instances and some preliminary computational results in Section 5. Section 6 includes the main computational study where we compare the proposed variants with the algorithms from the literature. Finally, we conclude the paper in Section 7.

Preliminaries
In this section, we provide the notation and definitions that are used throughout.
For p ≥ 2, R p is the p dimensional Euclidean space and · is the Euclidean norm on R p . The closed ball around a point a ∈ R p with radius r > 0 is denoted by B(a, r) := {y ∈ R p | y − a ≤ r}. Moreover, we use the following notation: e := (1, . . . , 1) T ∈ R p , e j ∈ R p is the unit vector with j th component being 1 and R p + := {y ∈ R p | y ≥ 0}. Let S be a subset of R p . The convex hull, conic hull, interior, closure, and boundary of S are denoted by conv S, co S, int S, cl S, and bd S respectively. A vector z ∈ R p \ {0} is a recession direction of S, if y + γz ∈ S for all γ ≥ 0, y ∈ S. The set of all recession directions of S is the recession cone of S and is denoted by rec S.
Let S ⊆ R p further be a convex set. A hyperplane given by {y ∈ R p | a T y = b} for some with y 1 , y 2 ∈ S and 0 < λ < 1 imply y 1 , y 2 ∈ F . A zero-dimensional face is an extreme point (or vertex ) and a one-dimensional face is an edge of S. A recession direction z ∈ R p \ {0} of convex set S is said to be an extreme direction of S if {v + rz ∈ R p | r ≥ 0} is a face for some extreme point v of S, see [21,Section 18].
This form is called a halfspace representation of S. If S has at least one vertex, then it can also be represented as S = conv V S + co conv D S , where V S ⊆ R p and D S ⊆ R p are the finite sets of vertices and extreme directions of S, respectively. For a polyhedral convex set S, we call two vertices adjacent if the line segment between them is an edge of S.
The problem of finding the set of all vertices and extreme directions of a polyhedral convex set S, given its halfspace representation is called the vertex enumeration problem. Throughout, we employ the bensolve tools for this purpose [17,18]. In addition to the vertices V and directions D of S, bensolve tools also yields the set of all adjacent vertices Let T ⊆ R p be another subset of R p . The Minkowski sum of S and T is S + T := {s + t ∈ R p | s ∈ S, t ∈ T }. Moreover, for α ∈ R, we have αS := {αs ∈ R p | s ∈ S} and we denote the set S + (−1) · T simply by S − T := {s − t ∈ R p | s ∈ S, t ∈ T }. Definition 2.1. Let S, T ⊆ R p . The Hausdorff distance between S and T is defined as where d(t, S) := inf s∈S t − s is the distance from point t to set S.
The following result provides a simple way of measuring the Hausdorff distance between a polyhedral set and a convex subset of it. The proof is omitted since it follows the same steps as the proof of [2, Lemma 5.3].
Lemma 2.2. Let S and T be convex sets in R p with rec S = rec T and T ⊆ S. If S is polyhedral convex with at least one vertex, then d H (S, T ) = max v∈V S d(v, T ) where V S is the set of all vertices of S.
Let C ⊆ R p be a nonempty closed convex cone. C is said to be pointed if C ∩ (−C) = {0} and non-trivial if C = {0}, C = R p . If C is pointed and non-trivial, then a partial order ≤ C is defined as follows: for all y 1 , y 2 ∈ R p , y 1 ≤ C y 2 holds if and only if y 2 − y 1 ∈ C.
The dual cone of C is defined as C + := {a ∈ R p | ∀x ∈ C : a T x ≥ 0} and it is a closed convex cone. If C ⊆ R p is a polyhedral cone, then C + is also polyhedral and can be written as C + = co conv {w 1 , . . . , w l }, where w 1 , . . . , w l ∈ R p are the extreme directions of C + for some l ≥ 0. In this case, y 1 ≤ C y 2 holds if and only if (w i ) T y 1 ≤ (w i ) T y 2 for all i ∈ {1, . . . , l}.
Let C ⊆ R p be a nonempty convex pointed cone. The following definitions are based on the partial order ≤ C and they are well-known concepts in convex vector optimization.
If C has nonempty interior and ({s}−int C)∩S = ∅, then s ∈ S is called a weakly C-minimal element of S. The set of all C-minimal (weakly C-minimal) elements of S is denoted by Min C S (wMin C S).

The problem
We consider a convex vector optimization problem given by where C is a closed convex pointed cone with nonempty interior, X ⊆ R n is a nonempty closed convex set and f : R n → R p is a continuous C-convex function given by f (x) := (f 1 (x), . . . , f p (x)) T . Throughout, we assume that C is polyhedral, in particular, its dual cone is given by C + := co conv {w 1 , . . . , w l }, and X is compact with nonempty interior. f (X ) := {f (x) ∈ R p | x ∈ X } is the image of X under f and the upper image for problem (P) is defined as P := cl (f (X ) + C), see Figure 1. P is a closed convex set and under the assumptions of the problem, P = f (X ) + C [2]. Moreover, the set of all weakly C-minimal points of P is bd P [7]. Different from a single objective optimization problem, there are various optimality and solution concepts for the problem (P). We will now introduce the solution concepts that will be used throughout. A feasible pointx ∈ X of (P) is called (weak) minimizer, if f (x) is (weakly) C-minimal element of f (X ). In the context of multiobjective optimization where the ordering cone is C = R p + , (weak) minimizers are referred to as (weak) efficient solutions. Another important concept for the multiobjective case that will be used through the rest of the paper is the ideal point, which is denoted by y I ∈ R p and defined as y I i := inf{f i (x) | x ∈ X } for i = 1, . . . , p.
Problem (P) is said to be bounded if P ⊆ {y} + C for some y ∈ R p . Since X is assumed to be compact, we know that problem (P) is bounded [16]. Below, we provide two solution concepts for bounded convex vector optimization problems, which are motivated from a set-optimization point of view [15]. The first solution concept depends on a fixed direction parameter c ∈ int C. It is introduced in [16] as a 'finite ǫ-solution'. Here we add the term Figure 1: Illustrations of sets f (X ), P (left) and solution concepts given by Definition 3.1 (middle) and Definition 3.2 (right). The dashed lines show bd (conv f (X ) + C) (middle and right) and solid blue lines show bd (conv f (X ) + C − ǫ{c}) (middle) and bd (conv f (X ) + C + B(0, ǫ)) (right). 'with respect to c' to emphasize this dependence. The second one is free of a direction parameter but depends on a norm in R p . . Let c ∈ int C be fixed. For ǫ > 0, a nonempty finite set X ⊆ X of (weak) minimizers is a finite (weak) ǫ-solution with respect to c if conv f (X ) + C − ǫ{c} ⊇ P.
Illustrations of an upper image and Definitions 3.1 and 3.2 can be seen in Figure 1. Note that both of these solution concepts yield an outer approximation to the upper image. On the other hand, conv f (X ) + C is an inner approximation to P. The Hausdorff distance between these inner and outer sets is bounded [6]. Moreover, as the feasible region X is compact, for any ǫ > 0, there exists a finite weak ǫ-solution (with respect to c) to problem (P), see [16,Proposition 4.3] and [2,Proposition 3.8].
There are different solution approaches to solve bounded CVOPs in the sense of Definitions 3.1 or 3.2. The main idea of these approaches is to generate (weak) minimizers for (P) in a structured way. One way of generating (weak) minimizers is to solve scalarization models. Below, we provide two well-known scalarization models together with some results regarding them.
For a weight parameter w ∈ R p , the weighted sum scalarization model is given by ). An optimal solution x ∈ X of (WS(w)) is a weak minimizer of (P) if w ∈ C + \ {0}. Conversely, for any weak minimizer x ∈ X , there exists w ∈ C + \ {0} such that x is an optimal solution of (WS(w)).
Pascoletti-Serafini scalarization [20] is given by The parameters v, d ∈ R p are referred to as the reference point and the direction, respectively. The Lagrange dual of the PS problem can be written as see [6,16]. Below, we provide some well-known results regarding (PS(v, d)) and (D-PS(v, d)).
The other implication follows similarly.
The following propositions show the existence of an optimal solution to problem (PS(v, d)) under the assumptions of problem (P) together with some additional conditions on the problem parameters. Note that Proposition 3.6 is already stated in [16,Proposition 4.4]. We still provide a proof of it as the one given in [16] has inaccuracies. 1 Proposition 3.6. Let v ∈ R p . If d ∈ int C, then there exist optimal solutions to (PS(v, d)) and (D-PS(v, d)). Moreover, the optimal values of the two problems coincide.
Proof. First, we show that there exists a feasible solution to (PS(v, d)). By the assumptions on (P), there existsx ∈ int X . Letz := max i∈{1,...,l} . . , l}. It is not difficult to see that (x,z) is a feasible solution for (PS(v, d)). Indeed, (x,z) is strictly feasible, hence a Slater point for (PS(v, d)). Hence, there exists an optimal solution w * for (D-PS(v, d)) and the optimal values of the two problems coincide.
For the existence of an optimal solution for (PS(v, d)), we first show that for any feasible solution (x, z) ∈ R n+1 , z is bounded below. As X is compact, (P) is bounded, that is, there exists a ∈ R p such that P ⊆ {a} + C. Then, by Remark 3.5, v + zd ∈ {a} + C holds for any feasible (x, z). This implies that for all w ∈ C + we have w T (v + zd − a) ≥ 0. Since w T d > 0 for all w ∈ C + \ {0}, it is true that Then, without changing the problem, we may add the constraint:z ≤ z ≤z to (PS(v, d)). As the feasible region of the equivalent form is compact, (PS(v, d)) has an optimal solution.
The next result states that using the primal-dual solution pair to problems (PS(v, d)) and (D-PS(v, d)), it is possible to find a supporting hyperplane to the upper image. It will play a critical role in the description of the solution algorithm which will be discussed in Section 4.  Let v, d ∈ R p and (x * , z * ), w * ∈ R p be the optimal solutions for (PS(v, d)) and its Lagrange dual, respectively. If the optimal objective function values of (PS(v, d)) and (D-PS(v, d)) coincides, then H :

The algorithm and variants
In this section, we will first provide the main framework of a CVOP algorithm that solves PS scalarizations iteratively and finds a finite weak ǫ-solution to problem (P) for a given ǫ > 0. Then, we will provide some variants based on selecting the parameters of the scalarizations.

The algorithm
We consider a CVOP algorithm that starts with an outer approximation to the upper image, iteratively updates it by solving PS scalarizations, and stops when the approximation is fine enough. More specifically, it starts by solving (WS(w i )) for all i ∈ {1, . . . , l} where w i 's are extreme directions of C + . Optimal solutions x i , i ∈ {1, . . . , l} of these problems form the initial set of weak minimizersX 0 . Then, the initial outer approximation P 0 of P is set as In the k th iteration, the vertices V k of the current outer approximation P k are considered and a vertex that was not used in the previous iterations, i.e., v ∈ V k \ V used is selected. We will later discuss different vertex selection methods, some of which return In this case, an upper bound for the Hausdorff distance between P k and P, namelyĥ = max v∈V k \V used z v can be computed. Ifĥ ≤ ǫ, the algorithm terminates by letting V k = ∅.
V k info = ∅ means that the vertex selection method does not store any information regarding the current vertices, see lines 6-12. If V k info = ∅ or if the algorithm is not terminated as explained above, then (PS(v, d)) is solved to find a weak minimizer x v , see Proposition 3.4. Note that if the direction parameter d is not fixed from the beginning of the algorithm, then it has to be computed before this step. The selected vertex v is added to V used and the corresponding weak minimizer x v is added to setX k (lines [13][14][15][16]. If the current vertex is close enough to the upper image, then the algorithm checks another vertex from V k \ V used . Otherwise, the current outer approximation is updated by intersecting it with the supporting halfspace H of P at f (x v ), see Proposition 3.8. The vertices of the updated outer approximation are computed by solving a vertex enumeration problem (lines [18][19][20]. The algorithm terminates if all the vertices of the current outer approximation are in ǫ distance to the upper image and returns a set of weak minimizersX .

Algorithm 1 Primal Approximation Algorithm for (P)
Compute the set of vertices V 0 of P 0 ; 5: Fix rules for selecting a vertex from V k and a direction parameter d; if V k info = ∅ or solve = 1 then 14: ). Let (x v , z v ) and w v be optimal solutions; 16: if z v > ǫ then 18: Compute the set of vertices V k+1 of P k+1 ;

21:
k ← k + 1; Later, we will discuss different rules for selecting the direction parameter and the vertices, respectively in Sections 4.2 and 4.3. The following proposition holds for any selection rule.

Proposition 4.2. When terminates, Algorithm 1 returns a finite weak ǫ-solution.
Proof. There exists optimal solutions x i to (WS(w i )) for all i ∈ {1, . . . , l} as X is compact and f is continuous. Moreover, holds since x i is an optimal solution for (WS(w i )) and inf c∈C (w i ) T c = 0. Then, P 0 ⊇ P. Moreover, P 0 has at least one vertex as C is pointed and (P) is bounded, see [21,Corollary 18.5.3]. Hence, we have x v is a weak minimizer. Hence, for any iteration k,X k is finite and contains only the weak minimizers (line 16). By Proposition 3.8, The algorithm terminates when V k \V used = ∅. Assume this is the case after K iterations. This suggests that each v ∈ V K is also an element of V used , that is, Similar to [2, Lemma 5.2], one can show that the recession cone of P k is C for any k. Hence, P K = conv V K + C holds true. Then, for any p ∈ P K , there exist ( . Hence, p ∈ conv f (X K ) + C + B(0, ǫ) and this implies (1) as p ∈ P K is arbitrary.
For the vertex selection rules that give V k info = ∅, assumeĥ = max v∈V K \V used z v ≤ ǫ holds for some K. Note that if there exists a vertex v in V K ∩ V used , then z v ≤ ǫ has to be satisfied by the structure of the algorithm. Hence, for the vertices V K of P K , we have z v ≤ ǫ. Similar steps for the previous case can be applied to show thatX K is a finite weak ǫ-solution of (P).

Direction selection rules
As explained in Section 1, algorithms from [16] and [7] follow the framework given by Algorithm 1. In [16], a fixed direction parameter d is used within (PS(v, d)) and in [7], d v is taken asp − v for a fixed pointp ∈ P. Before proceeding with the proposed direction selection rules, we use two examples (Examples 1 (for p = 2) and 2 provided in Section 5.1) to show that the selection of d in [16] andp in [7] affect the performance of these algorithms significantly. The results are summarized in Table 2, in which we report the number of scalarization models (SC) and the CPU time (T). 2 Table 2: Results for different d andp values fixed for algorithms in [16] and [7], respectively Example 1 (p = 2, ǫ = 0.005) Example 2 (ǫ = 0.05) d in [16]p in [7] d in [16] Motivated by these results, we propose two additional selection rules for the direction parameter to be used in PS scalarization.
Adjacent vertices-based approach (Adj). When the double description method is used to solve the vertex enumeration problem (line 20 of Algorithm 1), the sets of adjacent vertices for each vertex are also returned. We use this information to compute the direction parameter for a given vertex. In particular, for each vertex v, the normal direction of a hyperplane passing through the adjacent vertices of v is computed. The main motivation for this selection rule is that the positions of the adjacent vertices depend on the curvature of the upper image as they are formed by the supporting hyperplanes of the upper image at points that are possibly close to the current vertex. This geometric intuition is illustrated in Figure 3 for a two-dimensional setting.
Note that the adjacency list of a vertex may also contain extreme directions. If for a vertex v, we have an adjacent extreme direction z, then it is known that the line segment {v + rz ∈ R p | r ≥ 0} is a face of P k . In these cases, we construct an artificial adjacent vertex by moving the current vertex along the adjacent extreme direction, that is, we take v + z. An example of this case can be seen in Figure 3.
In general, a vertex may have more adjacent vertices than required. In that case, we adj , see line 3 of Procedure 2. Note that it is possible to obtain −d as the solution of this system. On the other hand, by Proposition 3.6, if d ∈ int C we obtain a weak minimizer by solving (PS(v, d)). Hence, we first check if one of two candidate unit normal vectors, d or −d, is in int C. If this is not the case, we use a predetermined direction d ∈ int C where d = 1. In . . , c r are the extreme directions of C. Remark 4.3. For p = 2, C = R 2 + , it is easy to show that either d ∈ int C or −d ∈ int C is satisfied (see Figure 3 for an illustration). However, this may not be true in general for Procedure 2 ChooseDirection(v,P ) Through our computational tests, we obtain some counterexamples for p = 4. Even though we haven't encountered this situation for p = 3, it is still possible that ±d / ∈ int R 3 + . See Figure 4 for an illustration, in which vertex v = 0 ∈ R 3 has 4 adjacent vertices: Ideal point-based approach (IP). For this approach, we assume that the ordering cone is R p + , hence the ideal point y I is well defined. For a vertex v of the current outer approximation, we consider the vector v − y I . Note that for the initial iteration, we have v − y I = 0. In the subsequent iterations, we obtain v − y I ≥ 0 since P k ⊆ P 0 for any k throughout the algorithm. For i ∈ {1, . . . , p}, we define d i =

Vertex selection rules
We propose different vertex selection rules that can be used within Algorithm 1. As will be detailed in Section 4.4, there are algorithms from the literature using some vertex selection rules that require solving additional optimization models, see [6,14]. Our motivation is to propose vertex selection rules which are computationally less complicated yet have the potential to increase the efficiency of Algorithm 1.
Vertex selection with clusters (C). This vertex selection rule clusters the vertices of the current outer approximation and visits these clusters sequentially. The motivation is to ensure that the vertices from different regions of the outer approximation are selected in a balanced fashion. The first step is to fix centers for the clusters, which can be done in different ways. For this purpose, we first solve (PS(v, d)) for each vertex of the initial polyhedron P 0 . The corresponding supporting halfspaces are intersected with the current polyhedron to obtain P 1 . The same procedure is repeated for P 1 and the vertices of P 2 are selected as the centers of the clusters, see Procedure 3. It is possible to use the vertices of P 1 to have fewer clusters or continue in the same fashion and use the vertices of P k for some k > 2 to have more clusters. Note that SelectCenters(V 0 ,P 0 ) has to be called right before line 6 of Algorithm 1.
For the remaining iterations of Algorithm 1, each vertex of the current outer approximation is assigned to the cluster whose center is the closest with respect to the Euclidean distance. Whenever a vertex has to be selected, a vertex is chosen arbitrarily from the cluster in turn. If there is no unexplored vertex assigned to the current cluster, then the algorithm selects the next nonempty cluster. The pseudocode is given in Procedure 4.
Vertex selection with adjacency information (Adj). Recall that bensolve tools returns the adjacency information for the vertices of the outer approximation. We use this to detect "isolated" vertices of the current outer approximation. The motivation is to obtain uniformity among the vertices of the outer approximation and consequently, the C-minimal points found on the boundary of the upper image as long as the geometry allows. For Compute the set of vertices V k+1 of P k+1 ; 7: end for 8: C = ∅. 9: for v ∈ V 2 do 10: t ← t + 1 and current ← t + 1 (mod |C|); 8: end while 9: Pick an arbitrary v ∈ C current ; 10: t ← t + 1; 11: return v, t each vertex v, the procedure finds the minimum distance, say dist v , from v to its current neighbors. Then, it selects the vertex which has the maximum dist v , see Procedure 5.
Vertex selection using local upper bounds (UB). This selection rule can only be used for C = R p + . It is motivated by split algorithms, which are originally designed to solve multiobjective integer programming problems, see for instance [5,12,3]. The main idea is to find a set of local upper bounds for the nondominated points and use them to select a vertex.
For this method, in line 1 of Algorithm 1, we additionally fix an upper bound u = Me where M is a sufficiently large number such that {u} − R p + ⊇ f (X ) and we initialize the set of upper bounds as U = {(u, ∅)}. Here, ∅ means that the upper bound u is not defined based on any other (weakly) nondominated point found in the algorithm. Note that the initial upper bound u satisfies that y I ≤ u. Through the algorithm, for any vertex v of P k , it is guaranteed that there exists a local upper bound u such that v ≤ u. Indeed, among all the upper bounds satisfying v ≤ u, we fix the one that yields the minimum v − u value as the 'corresponding' local upper bound for v.
Let V v adj be the set of vertices adjacent to v; 3: Together with the vertex v to proceed, this method returns also the corresponding local upper bound u where (u, y) ∈ U. Let (x v , z v ) be an optimal solution for (PS(v, d)). Using For the vertex selection procedure, we first assign the corresponding upper bounds for each vertex (line 3 of Procedure 7). Then, we choose v * ∈ V k \ V used , which has the greatest distance to its corresponding local upper bound u * (lines 10 and 12 of Procedure 7). Note that v * − u * yields an upper bound for the current approximation error since the corresponding upper bound of v. 3 As it will be detailed in Section 4.4, the algorithms proposed in [14] and [6] also compute an upper bound for the approximation error during any iteration. Different from them, UB does not solve optimization models to find this approximation error.
Note that some upper bounds have components equal to M. This causes the distances between v and its corresponding upper bound to be significantly large. To tackle this, the modifications are made, see lines 4-9 of Procedure 7 for the details.
Remark 4.4. For Procedure 7 to work correctly, we modify Algorithm 1 slightly as follows: In line 1, we initialize V 0 info as an empty set and in line 19 we also update V k+1 info ← V k info . The same modification is also applied to the algorithms from the literature that will be explained in Section 4.4.

Algorithms from the literature
In this section, we briefly explain similar approaches from the literature, namely, Algorithm KTW from [14], Algorithm DLSW from [6], and an algorithm from [2]. 5 The algorithm from [2] does not totally fit into the framework of Algorithm 1 as it solves a different scalarization model in each iteration. Here, we consider a variant of Algorithm 1, which is motivated by the one from [2].
Algorithm KTW. The outer approximation algorithm provided in [14] is a special case of the general framework given by Algorithm 1, applied to a maximization problem. Here, we shortly explain it for problem (P). KTW assumes 0 ∈ P. In iteration k, it minimizes the distance (with respect to an oblique norm) between the current outer approximation and the upper image. To do that, for every vertex v ∈ V k \ V used , it solves: Assuming (x v , λ v ) is an optimal solution to (2), v * ∈ arg min v∈V k \V used λ v is selected in is an optimal solution to (2) if and only if (x v , 1 − λ) is an optimal solution to (PS(v, −v)). Note that KTW is similar to the algorithm proposed in [7] in the sense thatp = 0 is fixed. However, different from [7], the selection of the vertices in each iteration is not arbitrary in KTW.
For the examples from Section 5.1, it is not guaranteed that 0 ∈ P. To overcome this, for the computational tests that will be presented in Section 6, we takep as in (6). Then, we modify the model given in (2) as: This simply corresponds to shifting the upper image as well as the vertices of the outer approximation byp ∈ P. To have a more efficient implementation to select the vertices, we do not solve (7) for all v ∈ V k \ V used . Instead, we check if (7) is solved in previous iterations.
Algorithm DLSW. Dörfler et al. [6] propose a vertex selection rule that uses the inner approximation obtained in each iteration. Accordingly, for each vertex v of the current outer approximation P k , the following problem, which measures the distance from v to the current inner approximation I k := cl conv f (X k ) + C, is solved: Then, vertex v * which is the farthest away from the inner approximation is selected, that is, v * ∈ arg max v∈V k \V used y v − v , where y v is an optimal solution of (8). Moreover, d An upper bound for the Hausdorff distance between the current outer approximation and the upper image is found as y v * − v * .
In [6], also an improved version of this algorithm is presented. Using the information from the previous iterations, it is possible to skip solving (8) for some of the vertices. In particular, for x k being the weak minimizer found in iteration k, if (8) is solved for some v in the previous iterations and (y v − v) T (f (x k ) − y v ) ≥ 0 holds true, then the solution of (8) for v in iterations k and k − 1 are the same by [6,Theorem 4.4]. Hence there is no need to solve (8) for v. See [6] for the details.
Algorithm AUU. Ararat et al. [2] propose an outer approximation algorithm to solve CVOPs. Even though this algorithm is similar to Algorithm 1, it is different since instead of (PS(v, d)), it solves the following scalarization: Note that (9) does not require a direction parameter. Instead, it computes the distance z v from v to the upper image, where (x v , z v ) is an optimal solution. If one solves (9) for all vertices of the current outer approximation P k , then it is possible to compute the exact Hausdorff distance between P k and P as max v∈V k z v , where V k is the set of vertices of P k . See [2] for the details.
As the main motivation of this study is to compare the effect of parameter selection rules for the PS scalarization, we consider a variant (namely, Algorithm AUU) of Algorithm 1 which is motivated by the algorithm from [2]. For AUU, we solve (9) only for parameter selection for the PS scalarizations and we still execute line 14 of Algorithm 1. The direction parameter is set as d = z v z v . Note that this is the direction that would yield the minimum optimal objective function value for (PS(v, d)).
In [2], the vertex selection is arbitrary, however, a special vertex selection rule is discussed in [2, Corollary 7.4 and Remark 7.5] to obtain a convergence result. Here, for algorithm AUU, we apply this vertex selection rule, that is, we solve (9) for each vertex of P k and select the one that is the farthest away from the upper image. This allows the algorithm to provide the approximation error in each iteration rather than at termination. Indeed, by Lemma 2.2, this approximation error is equal to the Hausdorff distance between P and the outer approximation. As in KTW, the solution of (9) for a vertex v is added to V k info , if (9) is not solved for v in the previous iterations.

Test examples and preliminary results
Before we proceed with the main computational study, we provide the numerical examples that will be used throughout. Moreover, to possibly decrease the number of variants to be tested for the main computational study, we provide some preliminary results.
For the computational tests provided here and in Section 6, the algorithms are implemented using MATLAB R2020b. The scalarizations are solved via CVX v2.2 [10,9] and SeDuMi 1.3.4 [23]. Vertex enumeration problems are solved by bensolve tools [18], which is a free solver providing commands for the calculus of convex polyhedra and polyhedral convex functions for Octave and MATLAB. The computer specification that is used throughout the computational study is Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz.

Test examples
Consider problem (P), where the ordering cone C is R p + ; and objective function f : R n → R p and feasible region X ⊆ R n are given as follows: Example 3. 7 Four instances of the following example are considered for a ∈ {5, 7, 10, 20}: Example 4. Three instances of the following example are considered for a ∈ {5, 7, 10}: For the computational tests, each component of A and P is generated according to independent normal distributions with mean 0, and variance 100; each component of b is generated according to a uniform distribution between 0 and 10.
Components are rounded to the nearest integer for computational simplicity. 6 This example is bounded in the sense that the upper image is included in y I + R 2 + where y I = 0 ∈ R 2 . In theory, the weighted sum scalarizations for the initialization step, namely, (inf x∈R+ x) and (inf x∈R+ We consider examples where the ordering cone is the positive orthant, mainly because some of the algorithms considered in our computational tests are designed to solve MOPs only. However, Algorithm 1 can solve vector optimization problems with ordering cones different from the positive orthant. To illustrate this feature, we consider Example 1 for p ∈ {2, 3} with different ordering cones taken from [2]. The details and results are provided in the appendix.

A preliminary computational study
As a preliminary analysis, we compare the proposed direction selection methods with the ones from the literature.
For the fixed point approach (FP) from [7],p ∈ P is set aŝ where x i is the optimal solution of (WS(w i )) for i ∈ {1, . . . , p} found in the initialization step and v 0 is a vertex of the initial polyhedron P 0 . Note that as the ordering cone is R p + , we have v 0 = y I . Moreover, it is not difficult to see thatp ∈ int P for all the test examples.
The direction parameter for a vertex v is set as d =p −v p−v . By Proposition 3.7, there are optimal solutions to (PS(v, d)) and (D-PS(v, d)) under this direction selection rule. However, for the computational tests, we ensure d ∈ int R p + for computational simplicity and to reduce the numerical issues. Accordingly, if d / ∈ int R p + then we use the predetermined direction d as in (Adj), which gives d = e e since C = R p + . For the fixed direction approach (FD) from [16], d is also set to e e . For our preliminary analysis, we employ different "simple" vertex selection rules which are summarized below.
• Closest to the Ideal Point (Id): Assume that the ordering cone is R p + and the ideal point y I is well-defined. Choose v ∈ V k \ V used such that v − y I is minimum.
• Farthest to an Inner Point (In): Choose v ∈ V k \ V used which yields the greatest distance p − v , wherep is computed as in (6).
• Random Choice (R): Randomly choose a vertex among V k \V used , using the discrete uniform distribution. For the numerical examples, we run the same example five times for this vertex selection rule.
For the test examples in Section 5.1, we run Algorithm 1 for direction selection rules FD, Adj, IP, FP and vertex selection rules Id, R, In. For the total of 12 variants formed by these parameter selection rules, we report the total number of scalarization models solved (SC) as well as the total CPU time (T) required, see Table 3. The results for Example 5 display the averages over 20 randomly generated problem instances.
We observe that while no direction selection method consistently outperforms the others for the same vertex selection method, FD and FP are outperformed by either Adj or IP in almost all cases except two: In Example 1 with p = 4 and vertex selection rule R, FP has the smallest SC; and in Example 4 with a = 7, FD has the smallest SC and T. In Examples 2-4, FP consistently performs worse than the other direction selection rules in terms of SC and T. On the other hand, in more than half of all the cases (21 out of 36) Adj performs the best among other direction selection rules, whereas IP performs the best in around one-third of all the cases (13/14 out of 36). Based on these results, we fix Adj as the direction selection rule for the main computational analysis.
When we compare the vertex selection rules for a given direction selection rule, we observe that R performs the best in almost all instances. The reason that Id and In decrease the algorithm's performance may be the result of these methods' tendency to choose vertices that are close to each other, especially for certain problem instances. However, R also has a disadvantage because it may not always provide consistent results.

Main computational results
In this section, we present a computational study on the examples listed in Section 5.1. We compare the proposed variants with the algorithms discussed in Section 4.4. We use three stopping criteria: the approximation error, CPU time, and cardinality of the solution set.
Recall that depending on our preliminary analysis from Section 5.2, we fix the direction selection rule as Adj. On the other hand, we consider all vertex selection rules, namely C, Adj, and UB as introduced in Section 4.3. In addition, we consider the random vertex selection, namely R, because of its good performance in our preliminary analysis in Section 5.2. Since the direction selection rule is common, we simply call the proposed variants C, Adj, UB, and R. For R, we solve each problem five times and report the average values.
Before presenting our computational study, we discuss two proximity measures. The first one is the approximation error, which is the realized Hausdorff distance between the upper image and the final outer approximation, i.e., d H (P K , P). Note that AUU returns the exact d H (P K , P) by its structure. For all the other algorithms, we compute the correct Hausdorff distance using the vertices v ∈ V K , after the termination.
The second measure is the hypervolume gap, which is originally designed in the context of multiobjective combinatorial optimization, see for instance [24], and recently used for convex vector optimization [1]. Here, we use it similarly to [1]. Let Q ⊆ R p be such that Q ⊇ f (X ) and V K be the vertices of the outer approximation at the termination of the algorithm. Let Λ be the Lebesgue measure on R p and Λ in := Λ((conv f (X K )+C)∩Q), Λ out := Λ((conv V K + C) ∩ Q). Then, we compute the hypervolume gap as HG = Λ out − Λ in .
Noting that the ordering cone is R p + for all the examples from Section 5.1, the set Q is taken as {u} − C, where upper bound u is set such that u i := max x∈X f i (x) for i = 1, . . . , p.
For the examples that u i cannot be computed (e.g., Example 2), the upper bound is set as where s is the number of the variants of Algorithm 1 that we solve the problem andX i is the solution set returned by the i th variant. For Example 2, the solutions found by the weighted-sum scalarizations at initialization are eliminated fromX i when computing u i . This is because these solutions are excessively large in one component and this causes bensolve tools not to perform the vertex enumeration as a result of numerical issues.
In the tables throughout, we report the total number of models solved to choose a vertex (VS), the number of scalarizations solved throughout the algorithm (SC), the cardinality of the solution set (Card), the CPU time (T), the realized approximation error (Err), and the hypervolume gap (HG). Note that VS is positive only for the algorithms in Section 4.4 as the others do not solve models to select vertices. The realized approximation error is d H (P K , P). For Example 5 with p = 3, we generate 20 instances as explained in Section 5.1 and report the averages of the results; whereas, for p = 4, we generate 3 instances and report the results of each instance separately. We use the convhulln() function in MATLAB to compute HG values. The ones that could not be calculated are given as (-) in the tables.

Computational results based on the approximation error
In this section, we solve the examples with a predetermined approximation error ǫ > 0, that is, when the algorithms terminate it is guaranteed that the Hausdorff distance between the outer approximation and P is less than ǫ. Tables 4 and 5 show the results for Examples 1-3 and 4-5, respectively. Since the realized errors are close to each other for all the algorithms, the tables do not show Err values. Moreover, since the number of scalarizations is the same as the cardinality of the solution set, Card is not reported separately.
For Example 4, we do not report the results for Algorithms AUU, DLSW, and KTW since bensolve tools was unable to perform the vertex enumeration in some iterations. For Example 5, DLSW results are not reported because of solver failures. Moreover, as the upper images are polyhedral for Example 5, the algorithms have the potential to return the exact upper image. We take ǫ sufficiently small (10 −8 ), hence the proximity measures are negligibly small and are not reported for these problems.  Tables 4 and 5 are analyzed, UB, R, C and Adj are faster than AUU, DLSW and KTW in all examples even though they solve more scalarization problems. Note that in Examples 1-3, AUU consistently solves the least number of scalarizations whereas, in Example 5, KTW yields the smallest SC values. Moreover, AUU has shorter CPU times  Table 4, we observe that AUU yields higher HG values than the others in most cases. The best HG values are obtained by DLSW or KTW, except in Example 1, p = 4, in which UB yields the smallest HG. Table 4 does not show the results for KTW for Example 2. Indeed, KTW can not solve this example, whose feasible region is not compact and upper image has an asymptotic behavior. Because of this structure of the upper image, depending on the position ofp, KTW finds larger distances through the iterations even though the approximations get finer. In order to illustrate this behavior, we pickp as (1.2, 1.2) T and perform a few iterations of KTW, see Figure 5. Note thatp is selected as a point that is close to the boundary of P, for illustrative purposes. This behavior of the algorithm may be observed for anyp, but generally after a certain number of iterations. In our computations, we did not observe this behavior illustrated in Figure 5 mainly because we fixp large enough. However, sincep is large, it caused numerical issues.

Computational results under limited runtime
In this set of experiments, we run Examples 1-5 under a runtime limit and report the proximity measures that the algorithms return in Table 6. For each example, around half of the shortest CPU time obtained in Section 6.1 is set as the time limit. For the variants UB, R, C and Adj, two different cardinality values are provided. The original value is greater (which is also equal to SC), whereas the smaller value is found by using Remark 4.1.
In Example 5, for p = 3, we solve the same randomly generated instances that we solved in Section 6.1. Since the CPU times for these instances differ significantly, the most timeconsuming 8 instances are selected, and the same runtime limit of 50 seconds is fixed for Figure 5: The outer approximations obtained after the first, second and third iterations of the KTW algorithm for Example 2 are shown, wherep = (1.2, 1.2) T . In each figure, the selected vertex for the corresponding iteration is indicated. The minimal point to be found after the scalarization is the intersection of bd P with the dashed line between the selected vertex andp. Even though the approximation gets finer in each iteration, the distances found by the KTW algorithm are 0.7619, 0.7633 and 0.7803, respectively.
them. The corresponding results in Table 6 are the averages over the 8 instances.
From Table 6, we observe that the solution sets returned by the variants (when Remark 4.1 is applied for UB, R, C and Adj) have similar cardinality values for two and three-dimensional examples. In Examples 2 and 3, C returned smaller approximation errors compared to the other variants, in general. In Example 1 with p = 3, proposed variants return smaller error values compared to the variants from the literature, but in Example 5 with p = 3, there is no clear conclusion.
For four-dimensional problems, AUU, DLSW, and KTW returned solution sets with smaller cardinality values than the others. In return, the approximation errors returned by AUU, DLSW, and KTW are larger than UB, C, and Adj for these problems. On the other hand, R behaves similarly to AUU, DLSW, and KTW. In most of the cases, HG values behave similarly to Err. However, in Example 5 for problems with p = 3 and (ins 3) of p = 4, the proposed variants return smaller HG, even though Err values are comparable among the algorithms.
An additional computational study is conducted on Example 3 to observe the behavior of the algorithms as the time limit changes. The variants are terminated after 50, 100 and 200 seconds. In this set of experiments, we also observe the effect of ǫ used within the variants UB, C, Adj and R even if the ǫ does not determine the stopping condition. Since we do not select the farthest away vertex to proceed with, ǫ has still an important role to decide if the algorithm would add a cut to the current outer approximation or not. To observe this effect, we run all the variants for the above-mentioned time limits both with ǫ = 0.005 and ǫ = 0. The results are provided in the appendix and the approximation errors obtained within different time limits for a = 5 are shown in Figure 6. 11 As expected, the performance of AUU, KTW and DLSW are not affected by the change in ǫ. On the other hand, an increase in ǫ affects the performance of UB, C and Adj positively. R is affected less compared to UB, C and Adj, while UB is affected significantly. Note that 10 ǫ is taken as in Tables 4 and 5 even though it is not the stopping criteria. 11 The results for other a values are similar, hence the corresponding figures are not included here. the two variants that give the worst approximation error in each setting (UB and R for ǫ = 0 and R and DLSW for ǫ = 0.005) are not plotted to increase visibility.

Computational results under fixed cardinality
In this section, we compare the performances of the algorithms when they run until finding a solution set with a predetermined cardinality. For each example, we give a cardinality limit so that the algorithms terminate before reaching the approximation error that is set for the experiments in Section 6.1. In Example 5 where p = 3, the same 8 instances are selected as in Section 6.2. For the variants UB, R, C, and Adj, we apply Remark 4.1. The results can be seen in Table 7.
In line with the observations from Section 6.1, we observe from Table 7 that UB, R, C, and Adj require less CPU time compared to AUU, DLSW, and KTW. Among the algorithms from Section 4.4, DLSW requires less CPU time compared to AUU and KTW, especially if the dimension of the objective space is high. When we compare the proposed variants, we see that Adj is slightly faster than the others in most examples.
The approximation errors obtained by AUU, DLSW, and KTW are smaller than or very close to those obtained by UB, R, C, and Adj, in general. This shows the tradeoff between the runtime and approximation error. When we compare the algorithms from Section 4.4, AUU and KTW yield better results in terms of the approximation error compared to DLSW, in general. On the other hand, the variants UB, R, C, and Adj are comparable among themselves. When we consider the HG values, we observe a similar behavior as Err, especially for problems with p ∈ {2, 3}. However, in Example 5 the proposed variants return smaller HG values even though the Err values are slightly worse than or comparable with

Conclusion
We consider bounded convex vector optimization problems and present a general framework for existing outer approximation algorithms from the literature. We observe that many algorithms iterate by solving Pascoletti-Serafini scalarizations or equivalent models and they mainly differ in selecting the parameters of the scalarization in each iteration. We propose additional methods to select the reference point (a vertex of the current outer approximation) and direction parameter of this scalarization. First, for a given vertex, we propose two methods to select the direction parameter. These methods are based on the position of the selected vertex as well as (1) the ideal point and (2) the adjacent vertices of the selected vertex, respectively. We compare the proposed direction selection rules together with the ones commonly used in the literature through preliminary computational tests and observe that using the positions of the adjacent vertices yields promising results. Hence, if the vertex enumeration problems are solved based on the double description method (e.g. by bensolve) so that the adjacency information is already obtained, then using Adj to select the direction parameter for each vertex is preferable.
We also propose some vertex selection rules which, different from the existing ones in the literature, do not require solving additional optimization problems. Instead, the proposed selection methods use the positions of all vertices (1) to form clusters and select accordingly, (2) to pick the most isolated vertex, or (3) to find the farthest away vertex from its corresponding local upper bound, respectively.
We implement the proposed vertex selection procedures together with the direction selection rule Adj and three relevant approaches from the literature. We provide an extensive computational study and observe that especially when the ordering cone is the positive orthant or a larger cone, the proposed variants perform better in terms of CPU time if the stopping condition is the approximation error or the cardinality of the solution set. On the other hand, the CPU times are comparable when the ordering cone is a strict subset of the positive orthant. 12 When the algorithms are run under a time limit, the proximity measures returned by the proposed variants are either better than or comparable to the ones returned by the algorithms from the literature. We also observe that the selection of ǫ affects the performance of the proposed variants even when ǫ is not used as the stopping criterion.
We conclude that if a reasonably good direction selection method (e.g. Adj) is used in the algorithm, then the vertex selection rule may have less impact on the performance of the algorithm. Hence, solving additional single objective optimization problems to select a vertex may not be necessary to improve the performance of the algorithm.
When we compare the overall performances of the proposed vertex selection methods, there is no clear conclusion. Indeed, random selection works very well in many cases, especially if the stopping condition is the approximation error. However, there are also some cases for which R works considerably worse than the proposed variants, see for instance the results for Examples 3 and 5 when the algorithms run until a predetermined runtime or cardinality limit (Tables 6, 7). In that sense, we conclude that R may not be as robust as the proposed variants, as expected. As a final remark, note that if the ordering cone is the positive orthant, then using UB has an advantage as it also returns an upper bound for the current approximation error when it is terminated after a predetermined runtime or cardinality limit. a computational study on Example 3 in which we compare three direction selection rules: d = e e (UB-FD), d = (u−v) u−v (UB-UB), and adjacent vertices approach (UB-Adj). We observe in Table 8 that UB-Adj solves less number of scalarizations while UB-UB clearly has a negative effect on T and SC. Based on this preliminary analysis, we confirm the usage of Adj as the direction selection rule throughout the computational studies in Section 6.
Let (x v , λ v ) be a solution to (7), Let z v be a solution to (9); 5:  Since algorithms KTW and UB are designed to solve multiobjective optimization problems, they are not included in this set of experiments. Table 9 shows the performances of AUU, DLSW, R, C, and Adj under different stopping criteria. Different from the main set of experiments, the bounding set Q considered for HG calculations is computed exactly as in [1, Section 6.1]. Also, these experiments are conducted using a different computer with specifications: Intel(R) Core(TM) i7-8565U CPU @ 1.80GHz. Table 9: Computational results for Example 1 for p = 2 with ordering cones C 1 and C 2 , and for p = 3 with ordering cones C 3 and C 4 . The first set (fixed ǫ) shows the results when the algorithms run until returning a finite weak ǫ-solution for provided ǫ values; the second set (fixed T) shows the results when the algorithms run under runtime limit T, and the third set shows the results when the algorithms run until finding a solution set with predetermined cardinality (Card).  Table 9, we observe that if the algorithms run until they return a finite weak ǫ-solution (fixed ǫ), then similar to the observations from Section 6, the proposed variants are faster than the algorithms from the literature for the ordering cones larger than the positive orthant (C 2 , C 3 ). However, there is no clear conclusion for the ordering cones smaller than the positive orthant. Under fixed runtime (fixed T), we observe that AUU and DLSW return solution sets with smaller cardinality. The algorithms are comparable in terms of the Err that they return. Moreover, AUU consistently yields the worst HG value and R, C and Adj return slightly better HG values than DLSW.
Finally, when the algorithms run until they find a predetermined number of solutions (fixed Card), the proposed variants return smaller HG values compared to AUU and DLSW. However, T and Err values are comparable for all algorithms.
D Results for different runtime limits for Example 3