Benders Decomposition for Bi-objective Linear Programs

In this paper, we develop a new decomposition technique for solving bi-objective linear programming problems. The proposed methodology combines the bi-objective simplex algorithm with Benders decomposition and can be used to obtain a complete set of extreme efficient solutions, and the corresponding set of extreme non-dominated points, for a bi-objective linear program. Using a Benders-like reformulation, the decomposition approach decouples the problem into a bi-objective master problem and a bi-objective subproblem, each of which is solved using the bi-objective parametric simplex algorithm. The master problem provides candidate extreme efficient solutions that the subproblem assesses for feasibility and optimality. As in standard Benders decomposition, optimality and feasibility cuts are generated by the subproblem and guide the master problem solve. This paper discusses bi-objective Benders decomposition from a theoretical perspective, proves the correctness of the proposed reformulation and addresses the need for so-called weighted optimality cuts. Furthermore, we present an algorithm to solve the reformulation and discuss its performance for three types of bi-objective optimisation problems.


Introduction
A wide variety of problems in different fields such as engineering and industry involve several conflicting objectives that should be taken into account simultaneously in optimisation problems.
It is generally impossible to find an optimal solution of such a multi-objective problem that optimises all the objectives concurrently.Many real-world decision making problems not only have multiple objectives, they are also modelled by complex mathematical programs that are trucks and drones to monitor forests for bushfires).Some researchers have integrated the concept of column generation into the steps of the biobjective parametric simplex method.Raith et al. (2012) describe how the entering variable selection in the bi-objective parametric simplex method can be reformulated as a column generation subproblem for BLPs.Moradi et al. (2015) extend this research to perform Dantzig-Wolfe decomposition by considering multiple subproblems, one for each commodity of a bi-objective multi-commodity network flow problem.
In this paper, we provide an introduction to bi-objective linear optimisation problems and the Benders decomposition technique in Sections 2.1 and 2.2.We present an approach that performs Benders decomposition for BLPs within the iterations of the bi-objective parametric simplex algorithm.Given a BLP, we first show that we can always find a bi-objective Benders reformulation consisting of only a subset of variables, two new variables to capture the first and second objective values of the remaining problem and cuts to ensure feasibility of the master variables as well as optimal objective values.This reformulation has a complete set of efficient solutions that is equivalent to that of the original BLP, see Section 2.3.Benders master and subproblems are formulated and an algorithm is developed that integrates the concept of Benders decomposition within the steps of the bi-objective parametric simplex method in Section 3.
Section 4 summarises the results of an extensive computational study.Section 5 concludes the paper and identifies promising directions for future work.This paper makes the following contributions to the literature.Firstly, we show how a full Benders reformulation can be derived in the bi-objective case, i.e. we introduce a bi-objective Benders reformulation of a non-decomposed BLP consisting of both feasibility cuts and what we term weighted optimality cuts.Secondly, we devise a bi-objective Benders decomposition approach and prove the correctness of the reformulation by showing that it can be used to obtain a complete set of efficient extreme solutions to the original problem.Finally, we present an algorithm, which we term the Bi-objective Benders Simplex Algorithm (BBSA), for solving the bi-objective Benders reformulation.This algorithm efficiently combines the bi-objective parametric simplex algorithm with Benders decomposition and enables solving the problem in a single sweep (with some backtracking).The BBSA is applied to solve three different types of BLP and extensive numerical results are discussed, including a comparison with a dichotomic approach that solves weighted single-objective Benders problems.

Benders Decomposition for Bi-objective Linear Optimisation Problems
In this section, we consider Benders decomposition in a bi-objective setting.To provide the necessary background, we begin in Section 2.1 by introducing a generic BLP and some relevant notation and terminology.In Section 2.2, we review how Benders decomposition can be used to solve a single-objective Linear Programme (LP).Section 2.3 then extends the single-objective methodology to the bi-objective case.In particular, we introduce the bi-objective Benders reformulation of a BLP supported by an illustrative example.
Figure 1: Feasible set Z in objective space, set of non-dominated points Zn (black dotted) and non-dominated extreme points z 1 , z 2 , z 3 , z 4 .

Bi-objective Linear Programs
A BLP in standard form is given by min (ĉ 1 ) ⊤ x where x ∈ R n and ĉ1 , ĉ2 , Â, b have appropriate dimensions.In the following we define some terms related to solutions of BLPs.
Definition 1.For a BLP we call R n the decision space and R 2 the objective space.The set of feasible solutions in decision space x ∈ R n that satisfy the constraints of BLP is called the feasible set in decision space X ⊆ R n , whereas its image is called the feasible set in objective space Z ⊆ R 2 .For (1) the feasible set is X = {x ∈ R n | Âx = b, x ≧ 0} and its image in objective space is Z = z(x) ∈ R 2 |z(x) = (z 1 (x), z 2 (x)) ⊤ and x ∈ X .
Definition 2. When comparing objective vectors z, z ∈ R 2 the following notation is used.We define symbols <, ≤, and ≦ as follows: z < z if z i < zi for i = 1, 2; z ≦ z if z i ≦ zi for i = 1, 2; z ≤ z if z i ≦ zi for i = 1, 2 and z ̸ = z.When z ≤ z, z dominates z, and when z < z, z strictly dominates z.
Definition 3. When solving a bi-objective optimisation problem we seek to identify the set of efficient solutions X e that consists of all solutions x ∈ X whose objective vector z(x) is not dominated by that of another solution x ∈ X .The corresponding set Z n in objective space is called the set of non-dominated points Z n = {z(x) ∈ R 2 |x ∈ X e }.Similarly, a weakly efficient solution is defined as a solution x ∈ X whose objective vector z(x) is not strictly dominated by that of another solution x ∈ X .We aim to identify a complete set of efficient (extreme) solutions that contains at least one efficient (extreme) solution per non-dominated point z ∈ Z n .
The feasible set in decision space of an LP (and a BLP) is a polyhedron, and so is the feasible set Z of BLP in objective space, see also Figure 1.The set of non-dominated points Z n is shown as a dotted line in the figure, and it lies on the lower-left boundary of Z.
It is known that a feasible solution x ∈ X of BLP ( 1) is efficient if and only if there exists a weight λ ∈ (0, 1) such that x is an optimal solution of λ-(1), Isermann (1974).As illustrated in Figure 1, all non-dominated points lie on facets of the polyhedron in objective space, where any non-dominated point that is not an extreme point corresponds to an optimal solution of λ-(1) for exactly one value of λ.Some non-dominated points are optimal for a range of λvalues and lie on the intersection of adjacent facets, such as points z 1 , z 2 , z 3 , z 4 in Figure 1.
These non-dominated points are termed non-dominated extreme points and they correspond to efficient extreme solutions.Note that, for a problem with bounded set of efficient solutions (Remark 1), all non-dominated points of a BLP can be obtained as convex combinations of neighbouring non-dominated extreme points, and therefore a complete set of efficient solutions (and corresponding non-dominated points) can be obtained by taking all convex combinations of neighbouring efficient extreme solutions (and corresponding non-dominated extreme points).
The non-dominated extreme points are z 1 , z 2 , z 3 , z 4 in Figure 1.
To solve a BLP by weighted-sum scalarisation one must iteratively solve single-objective weightedsum problems to obtain a complete set of efficient extreme solutions, and methods such as the dichotomic approach (e.g.Aneja and Nair, 1979;Cohon et al., 1979) provide a streamlined approach for this.The bi-objective parametric simplex algorithm is an alternative approach that does not iteratively solve single-objective weighted-sum problems.This method determines a complete set of efficient extreme solutions X e = {x 1 , x 2 , . . .x q } with increasing first objective value, i.e. z 1 (x 1 ) < z 1 (x 2 ) < . . .< z 1 (x q ).Each efficient extreme solution x i is an optimal solution of λ-BLP for all λ ∈ [a i+1 , a i ] with a i ≧ a i+1 and a 1 = 1, a q+1 = 0 for i = 1, . . ., q.
We provide an overview of the bi-objective parametric simplex algorithm and further details can be found in Ehrgott (2005).Through a sequence of basis updates, the algorithm iteratively moves from one non-dominated point to the next, which corresponds to proceeding from left to right in Figure 1.The algorithm starts by finding an efficient extreme solution x * that is optimal with respect to the first objective function.Associated with this solution is a set of basic variable indices and a set of non-basic variable indices.As per the standard simplex algorithm, the bi-objective version iteratively performs basis exchanges by identifying entering and leaving variable pairs.In a bi-objective setting, the entering variable is identified to be a non-basic variable with maximum improvement of the second objective function over the deterioration of the first objective function.This ratio can be computed using the reduced costs, c k , for objective functions k = 1, 2. If an entering variable exists, then the leaving variable can be determined (if it exists) as per the simplex algorithm for single-objective LP.This iterative procedure of identifying an entering variable terminates when no entering variable can be found (i.e. the basis is optimal with respect to the second objective function).
Remark 1.We assume X ̸ = ∅.In general the set of efficient solutions of a BLP is made up of efficient extreme solutions, and / or unbounded efficient edges, or rays, representing directions in which the set of non-dominated points in objective space is unbounded (Rudloff et al., 2017).
We also assume that the set of efficient solutions X e of all considered bi-objective optimisation problems is bounded (no unbounded efficient edges exist).This means that the weighted-sum scalarisation of the BLP is also bounded, i.e. that an optimal solution of λ-(1) exists for λ ∈ [0, 1].
It should be noted that the parametric simplex algorithm (and our proposed algorithm) can be adapted to the unbounded case for BLP following Rudloff et al. (2017).

Single-Objective Benders Decomposition
Problem size is critical when solving optimisation problems.In many practical applications, the number of decision variables or constraints is prohibitively large.Benders decomposition is a technique that can be used to deal with large-scale LPs.As this paper extends this approach to a bi-objective setting, we provide a comprehensive overview of the approach.Further details can be found in e.g., Benders (1962); Taskin (2010).Consider a single-objective LP of the form where again x ∈ R n , y ∈ R q and vectors c, f, b, d and matrices A, B, D have appropriate dimensions.One can observe that (3) reduces to a problem in terms of the x variables once the the y variables have been fixed to certain values, and can therefore be restated as (4) where θ(ȳ) is the optimal objective value of (5), which must be considered for all ȳ ∈ Y, where If ( 5) is bounded, then the value of θ(ȳ) can also be obtained by optimising its dual formulation.
Assuming that π denotes the vector of dual variables associated with the constraints of (5), then the corresponding dual formulation can be stated as: An advantage of solving the dual formulation is that its feasible region is independent of the values of the y variables.The feasible region D of ( 6) can be stated in terms of its set of extreme points D p and extreme rays D r , i.e.
Therefore (6) can be restated to give: where θ ∈ R is the only decision variable.The first set of constraints, termed feasibility cuts, ensures that any ȳ ∈ Y that results in an unbounded direction for ( 6) is infeasible.The second set of constraints, termed optimality cuts, ensures that the θ variable correctly captures the optimal objective value of the x variables.Using (7), the following Benders reformulation of the original problem (3) can be obtained.
Practical applications of Benders decomposition iteratively identify and add violated feasibility and optimality cuts rather than enumerating the sets D r and D p a priori, which may be prohibitively large.The master problem is a relaxation of (8) that includes only subsets of feasibility and optimality cuts.An iterative approach finds candidate solutions (y, θ) that are then subsequently checked for infeasibility/optimality by the subproblem (6).If ( 6) is unbounded, a feasibility cut is added to the master problem, while, if (6) is feasible, a check is made to see whether or not y is optimal.An optimality cut can be generated and added to the master problem to remove the provably suboptimal solution y if θ(y) > θ.Lower and upper bounds on the optimal objective value of the master problem are obtained at each iteration, and the problem is solved once the bounds are considered sufficiently close in value.

Benders Reformulation of BLP
We now focus on developing a bi-objective variant of Benders decomposition and consider BLPs of the form where x ∈ R n , y ∈ R q and vectors c 1 , c 2 , f 1 , f 2 , b, d and matrices A, B, D have appropriate dimensions.Note that BLP (1) was originally formulated in its standard form, whereas we use the more convenient inequality form (9) here.To apply Benders decomposition to this problem in which we leave the y variables in the master problem and place the x variables in the subproblem, we proceed as in Section 2.2 and reformulate (9) as a Bi-objective Master Problem (BM) in (10) where θ 1 (ȳ) and θ 2 (ȳ) bound the objective values of the Bi-objective Subproblem (BS) in ( 11), which must be considered for all ȳ ∈ Y: It should be noted that (11) will generally have infinitely many efficient solutions in the decision space and non-dominated points in the objective space.Variables θ 1 and θ 2 are used to capture the contribution of variables x to the two individual objective functions in (10), enabling the correct representation of non-dominated points in objective space.
How to handle a bi-objective subproblem (11) in the context of Benders decomposition is an interesting question to ask.We first consider the impact of using two separate subproblems, one for each objective, and investigate the conditions under which such an approach is sufficient for obtaining a complete set of efficient extreme solutions to BLP (9).Both subproblems are of the form (5) with the exception that the first (respectively, second) has the objective function (c 1 ) ⊤ x (respectively, (c 2 ) ⊤ x).Continuing as in the single-objective case, we identify the set of extreme points for the dual formulation of each subproblem.These are denoted D 1 p and D 2 p , respectively.Note that the set of extreme rays to each of the two dual formulations is the same as the set of rays D r is independent of the subproblem objective in (5).This leads to the following Benders reformulation: If one then proceeded as in the single-objective case (8), dynamically separating violated optimality and feasibility cuts using two subproblems instead of one, one might expect that (12) correctly describes a complete set of efficient extreme solutions of the original problem (9).Our first observation is that (12), which only relies on single-objective subproblems, is not able to fully describe all efficient solutions, as the following example illustrates.
We apply a Benders decomposition in which the y variables are placed in the master problem, and the x variables are placed in the two single-objective subproblems.The feasible regions of the resulting dual subproblem formulations are illustrated in Figure 2; the blue (respectively, green) area in the left (respectively, right) plot corresponds to the first (respectively, second) objective.
There are also two extreme rays D r = {(1, 0), (0, 1)} for each space.The Benders reformulation can therefore be stated in the form of (12) as Reformulation (14) only has three non-dominated extreme points {z 1 , z 2 , z 5 }.It correctly identifies points z 1 and z 2 as shown in Figure 3 but incorrectly identifies point z 5 = − 2 3 , − 18 5 .This point is defined by the first objective value for point z 3 and the second objective value for point non-dominated points of (13) non-dominated points of ( 14) z 4 .Points z 3 and z 4 have identical master variable values (y 1 , y 2 , y 3 ) = 0, 1 3 , 0 but different subproblem variable values (x 1 , x 2 ).
This example demonstrates that it is insufficient to only add optimality cuts corresponding to extreme dual solutions to the single-objective subproblems.To understand why, we analyse the efficient extreme solution associated with point z 4 .In this solution, the subproblem variables take the values (x 1 , x 2 ) = 1 5 , 3 5 .The corresponding dual solutions are π 1 1 , π 1 2 = 1 6 , − 13 15 for the first subproblem and π 2 1 , π 2 2 = 1 3 , 4 15 for the second, respectively.We note that the first dual solution is infeasible.The explanation for this is that the solution (x 1 , x 2 ) = 1 5 , 3 5 is not optimal for the first objective.The second dual solution is feasible (and hence optimal for the second objective).To ensure θ 1 accurately captures the relevant subproblem objective function value, the following cut could be defined: This is not present in (12) since it is derived from a dual solution that is not an extreme point of the dual feasible region for the subproblem associated with the first objective.Evaluating this constraint at the master solution (y 1 , y 2 , y 3 ) = 0, 1 3 , 0 enforces θ 1 ≥ 1 5 , the correct bound for the first objective function value in z 4 .However, the bound implied in (14) at the same master solution is θ 1 ≥ − 2 3 , (due to the single optimality cut for θ 1 ).This is the correct bound for the first objective function value in z 3 .Being more restrictive, the addition of (15) to (14) would render it impossible to generate the third non-dominated extreme point (z 1 , z 2 ) = − 2 3 , − 10 3 .From this we conclude that single-objective optimality cuts alone are insufficient.
A weighted cut, which weights the values of θ 1 and θ 2 in an optimality cut, would allow θ 1 and θ 2 to take different values for the same master solution (y 1 , y 2 , y 3 ) = 0, 1 3 , 0 .Such a cut enables the Benders reformulation to correctly identify the two efficient extreme solutions with (x 1 , x 2 ) = 1 5 , 3 5 and (x 1 , x 2 ) = 0, 2 3 , respectively.We demonstrate this by adding a weighted cut with λ = 4 17 .This weight is specifically chosen since it defines the weight for which the non-dominated extreme points z 3 and z 4 are optimal for λ-(13).The weighted cut is of the form: where (π 1 , π 2 ) = 5 17 , 0 is the optimal dual solution of the λ = 4 17 weighted subproblem that is induced by the master problem solution (y 1 , y 2 , y 3 ) = 0, 1 3 , 0 .Adding just this cut to reformulation (14) yields the correct set of efficient extreme solutions (and hence non-dominated extreme points), as shown in Figure 3.
Then there exists λ ∈ (0, 1) (Isermann, 1974) such that (x, ȳ) is an optimal solution of the weighted-sum problem (9) with objective function (c Applying this weighting factor λ to BM (10) with c λ , f λ as before, leads to a weighted (primal) subproblem ( 16), which we will also refer to as λ-(11): The feasible region of the dual of λ-( 11) can be described as ).The (single-objective) Benders reformulation (see also ( 7)) of the λweighted problem λ-( 9 Definition 4. Given λ ∈ [0, 1] and a corresponding π λ p ∈ D λ p .Then a weighted optimality cut takes the following form: We therefore obtain the following Benders reformulation of BLP (9), which includes weighted cuts of the form (17) with λ ∈ [0, 1].As such, the single-objective cuts corresponding to λ = 0 and λ = 1 are included in the set of weighted cuts.
(18) Formulation (18) represents the complete Benders BLP reformulation.While there is an infinite number of weighted cuts, we will show that only a finite set of them is required in (18), and that it is always possible to generate this finite Benders BLP reformulation.We use (18)-{ Dr , Λ, Dp } to refer to a variant of (18) that replaces D r by Dr and D λ p for λ ∈ [0, 1] by Dp (λ) for λ ∈ Λ, with Dr ⊆ D r and Dp (λ) ⊆ D λ p .In the following, we prove that, while the efficient solutions in decision space of ( 9) and ( 18) may be different, their respective sets of nondominated points in objective space are equal.Furthermore, the reformulation only requires a finite subset of cuts Dr and Dp (λ) for λ ∈ Λ.
Proof.We begin by observing that the constraint set of ( 18) is the union of the constraint sets of an infinite number of Benders reformulations of the single-objective λ k -( 9) problem.The same weighted-sum scalarisation can be applied to the two objective functions of ( 9) and ( 18) for some λ k ∈ [0, 1].From the theory on single-objective Benders decomposition we know that the optimal objective value of λ k -( 9) is equal to the optimal objective value of λ k -( 18).Consequently, for any λ k ∈ [0, 1] the following relationship between an optimal solution (x ′ , y ′ ) of λ k -( 9) and an optimal solution (θ ′′ 1 , θ ′′ 2 , y ′′ ) of λ k -(18) exists: where , where x ′′ is an optimal solution to min {(c λ k ) ⊤ x|Ax ≧ b − By ′′ , x ≧ 0}, i.e., the subproblem ( 16) that is induced by y ′′ .
Each efficient extreme solution (x i , y i ) ∈ X e for i = 1, 2, . . ., q is an optimal solution to λ k -( 9) This must be an optimal solution to λ k -( 18) for all λ ∈ [a i+1 , a i ].If this were not the case, then for some i = 1, 2, . . ., q, there must exist a λ j ∈ (a i+1 , a i ) and an optimal solution (θ j 1 , θ j 2 , y j ) to λ j -(18) for which However, this contradicts the observation in ( 19).If the solution (θ j 1 , θ j 2 , y j ) were an optimal solution to λ j -( 18), then there must exist a feasible solution (x j , y j ) to λ j -( 9) with the same objective value.This contradicts the assumption that (x i , y i ) is an optimal solution to λ k -( 9) for all λ k ∈ [a i+1 , a i ].We can conclude that (θ j 1 , θ j 2 , y j ) does not exist.Each solution (θ i 1 , θ i 2 , y i ) ∈ X r e for i = 1, 2, . . ., q is therefore an optimal solution to λ k -( 18) for all λ k ∈ [a i+1 , a i ] and has individual objective function values z i 1 and z i 2 .We also note that any y i , where i = 1, 2, . . ., q, can be used to provide an upper bound to λ k -(9), and hence λ k -(18), with i ̸ = k, since the feasibility of y i is independent of the objective function.Feasibility cuts, defined by D r , do not depend on the λ k value.
For a given λ k ∈ [0, 1] the constraints of (18) ensure that any feasible solution cannot have a better λ k -weighted objective value than the optimal objective value of the weighted single objective problem λ k -(9).Thus, any solution that is efficient in (9) will have a corresponding efficient solution in (18).This is also true for efficient extreme solutions.Thus, X r e is a set of efficient extreme solutions to (18).The set of non-dominated extreme points remains Z n with this reformulation.Since, each solution (θ i 1 , θ i 2 , y i ) is optimal for the range λ ∈ [a i+1 , a i ] we can remove all λ ∈ (a i+1 , a i ) for i = 1, 2, . . ., q without changing the set X r e .Therefore we can choose Λ = {a 1 , a 2 , . . ., a q }, Dr = D r , and Dp (λ) = D λ p for λ ∈ Λ to obtain (18)-{ Dr , Λ, Dp }.
We conclude this section by identifying conditions under which single-objective optimality cuts, defined by D 1 p and D 0 p , are sufficient to define the Benders reformulation of (9).
Remark 2. Consider a Benders reformulation of (9) of the form (18), with master variables y and subproblem variables x.If the optimal solution to the weighted subproblem (16) is the same for λ = 1 and λ = 0 for any master (10) solution ȳ, i.e, if the two subproblem objectives in (11) are not contradictory, then weighted cuts, obtained using any λ k ∈ (0, 1) are not necessary in the Benders reformulation.
Please refer to Appendix A for a proof of Remark 2.
Remark 3. Consider a Benders reformulation of (9) of the form (18), with master variables y and subproblem variables x.Assume that the subproblem variables are only present in one objective.If, w.l.o.g., this is the first objective, we have c 2 = 0 in (9).Therefore, the subproblem (11) reduces to a single-objective optimisation problem, variable θ 2 can be omitted from (18) and, again, only single-objective cuts for λ = 1 are required.

Algorithm for Bi-objective Benders Decomposition
Proposition 1 confirms that for a BLP (9), a bi-objective Benders reformulation (18)-{ Dr , Λ, Dp } with a finite number of weighted cuts exists.To obtain such a reformulation one could iteratively solve single-objective weighted-sum problems with Benders decomposition, for instance by applying a dichotomic approach (see also Section 2.1).This scalarisation approach would yield a single-objective Benders reformulation (8) of each weighted problem λ-(9) for some λ ∈ [0, 1].
The cuts from all scalarised reformulations with λ ∈ Λ are combined to obtain (18)-{ Dr , Λ, Dp } where the optimality cuts in (8) are converted to weighted cuts by replacing . Instead, we describe an algorithm in the following that integrates the generation of Benders feasibility and weighted optimality cuts into the parametric bi-objective simplex algorithm to solve BM (10) and BS (11), thereby not solving single-objective scalarised problems.
The two approaches are compared in Section 4.8.
In this section, the proposed BBSA is introduced.The BBSA works by iteratively identifying non-dominated points of a Benders master problem BM by applying bi-objective parametric simplex iterations, where the initial BM is a problem of the form (10), which is (18)-{∅, ∅, ∅} obtained from ( 18) by removing all feasibility and optimality cuts.At any iteration of the algorithm, this Benders reformulation is made up of a subset of feasibility cuts Dr and optimality cuts Dp (λ) for λ ∈ Λ that have been added up to the current iteration of the BBSA, and is the problem we refer to as BM (18)-{ Dr , Λ, Dp } in the context of the BBSA.BM (18)-{ Dr , Λ, Dp } is iteratively updated by adding all identified feasibility and optimality cuts.
We term an efficient solution ( θ1 , θ2 , ȳ) with objective vector z of BM ( 18 where π r ∈ D r and Dr is updated to include π r (i.e.Dr = Dr ∪ {π r }).Otherwise, the values of θ 1 and θ 2 in BM (18)-{ Dr , Λ, Dp } are modified, if necessary, by adding one or more weighted optimality cuts for π λ i p ∈ D p (λ i ), as explained later in this section, to obtain an updated BM (18)-{ Dr , Λ, Dp } (i.e.Λ and Dp are updated).The BBSA terminates when there are no remaining unexplored non-dominated points of BM, at which point the final (18)-{ Dr , Λ, Dp } is obtained.
The BBSA investigates at most one unexplored non-dominated point of BM (18)-{ Dr , Λ, Dp } in each iteration.It is therefore not necessary to fully solve BM (18)-{ Dr , Λ, Dp } with the bi-objective simplex algorithm in each iteration of the BBSA.Let us suppose that, at a given iteration, amongst all non-dominated points of BM (18)-{ Dr , Λ, Dp }, z last is the most recently explored (and therefore confirmed) point and that (θ last 1 , θ last 2 , y last ) is the corresponding efficient solution.In the subsequent iteration, the bi-objective simplex algorithm starts from the explored non-dominated point z last and moves to a neighbouring unexplored non-dominated point.
Algorithm 1 provides an overview of the BBSA algorithm.The input of Algorithm 1 is problem data, and its output is the set of non-dominated extreme points Z nx for problem (9) and (18)-{ Dr , Λ, Dp }, and corresponding efficient solutions X ex of (18)-{ Dr , Λ, Dp }.Sets Zn and Xe track non-dominated points and efficient solutions found while the algorithm runs (initialised in line 1).Weight set Λ and weight λ = 1 are also initialised in line 1.
The BBSA starts by finding the first non-dominated extreme point, which is optimal for the first objective function, by solving BM i.e. λ-(18)-{∅, ∅, ∅} for λ = 1 and without cuts, for the first objective function using the standard Benders decomposition technique, see line 2, function SolveSingleObjectiveBenders.Solution ( θ1 , θ2 , ȳ) and objective cost vector z are returned as well as the resulting sets of extreme rays and extreme points denoted Dr , Dp (1) that form the initial cuts for problem (18)-{ Dr , Λ, Dp }.
The BBSA iterates between exploring a solution in line 4, i.e. solving the subproblem BS (11) with the parametric simplex method (if it is feasible) and adding cuts (lines 5-13), storing solutions that are confirmed as efficient (lines 18-21) and partially solving the updated BM Algorithm 1: BBSA to solve bi-objective linear optimisation problems Input: ȳ and (implicitly) (lines 23-31).In the following we describe these steps in detail.
In line 4, a solution ȳ is explored by applying function Explore (Algorithm 2).Inputs to Algorithm 2 are ȳ and the problem parameters required to define BS (11).It returns either a feasibility cut D + r or weighted optimality cuts D + p , Λ + to be added to BM.At first the BLP BS ( 11) is solved in line 1 of Algorithm 2 where function SolveBLP applies the bi-objective parametric simplex (e.g.Ehrgott, 2005), returning the set of k efficient extreme points X BS = {x 1 , x 2 , . . ., x k } of BS as well as the corresponding non-dominated extreme points Z BS = {z 1 BS , z 2 BS , . . ., z k BS }.Set A = {a 1 , a 2 , . . ., a k , a k+1 } with a k+1 = 0 contains the weight limits such that for } is also returned and contains, for each extreme point x i ∈ X BS , vectors of dual variables (π i 1 , π i 2 ) associated with the constraints of BS (11), the values of which are determined by the respective objective function.Vector π r is an unbounded extreme ray (or a null vector if none exists).If the problem is infeasible (X BS = ∅), only the extreme ray is returned in D + r in line 4 and the others sets remain empty (lines 5 and 6).If, on the other hand, the problem is feasible, D + r remains empty (line 9), and the weighted duals are returned as D + p (line 10) and corresponding weights are Λ + (line 11).We note that it is assumed the ordering of elements in X BS , Z BS , Π, A is maintained throughout the algorithm.
In Algorithm 1 the cut sets are updated next.If the BS (11) was found to be infeasible in Explore(ȳ), Dr is updated by adding a feasibility cut for the identified extreme ray (line 6).Otherwise the BS (11) was found to be feasible and, for each weight a i ∈ Λ + , a weighted optimality cut (17) is added (lines 9-11), where λ = a i , and the weighted dual is Note that an additional optimality cut for θ 2 is generated in line 12 for λ = 0 in (17).Here, is the vector of duals corresponding to the last non-dominated BLP (9).Finally, Proposition 3 confirms that the algorithm terminates after a finite number of steps.
Lemma 1.A solution ( θ1 , θ2 , ȳ) ∈ X ex with image z ∈ Z nx identified by Algorithm 1 always corresponds to an efficient solution of (9), i.e. there exists an efficient solution of (9) with non-dominated image in objective space that is identical to z.
If BS ( 11) is feasible, it is fully solved in line 4 of Algorithm 1 (with details in Algorithm 2) with all weighted cuts added in lines 9-12.For each efficient extreme solution of BS, there is an interval of weights [a k+1 , a k ] for which the efficient solution is an optimal solution for λ-BS, In lines 9-12 a weighted cut is generated for each obtained efficient extreme solution of BS and corresponding λ ∈ {a 1 , a 2 , . . ., a k+1 }, with a 1 = 1 and a k+1 = 0.In particular, there cannot exist a weighted cut for λ ∈ [α, β] that would render solution ( θ1 , θ2 , ȳ) infeasible.
To see this choose any λ ∈ [α, β], which will be a convex combination λ = γa k+1 + (1 − γ)a k , γ ∈ [0, 1] of the interval limits of one of the weight intervals [a k+1 , a k ] identified when solving BS.In particular the following two weighted cuts were added for a k+1 and a k (amongst others): If solution ( θ1 , θ2 , ȳ) is feasible for weighted optimality cuts ( 21) and ( 22), it is also feasible for their convex combination as the following is a convex combination of ( 21) and ( 22): where π λ p is an optimal dual solution of the weighted subproblem λ-(11) (see also definition of a weighted optimality cut in ( 17)).Please refer also to the proof of Remark 2 in Appendix A, which uses a similar argument.
We can conclude that the solution ( θ1 , θ2 , ȳ) is an optimal solution to the corresponding prob- Proposition 2. Algorithm 1 terminates with a bi-objective master problem BM (18)-{ Dr , Λ, Dp } whose set of non-dominated extreme points Z nx is equal to that of the original BLP (9).
Proof.Observe that at any iteration of the BBSA, the master problem that is considered, BM (18)-{ Dr , Λ, Dp }, constitutes a relaxation of (18) with Dr ⊆ D r and ∪ λ∈Λ Dp (λ) ⊆ ∪ λ∈[0,1] D λ p .Being a relaxation, any feasible (θ 1 , θ 2 , y) of ( 18) is also a feasible solution of BM (18)-{ Dr , Λ, Dp }; however, the reverse is not necessarily true.We note that, in particular, the set of efficient solutions X e to (18) remains feasible for (18)-{ Dr , Λ, Dp } at every iteration.Given a solution (θ 1 , θ 2 , y) ∈ X e there does not exist a i.e., no violated, weighted optimality cut exists.The BBSA cannot identify any dual extreme point of λ-( 11), for a particular choice of λ ∈ [0, 1], that is not contained in ∪ λ∈[0,1] D λ p , since the latter is an enumerated set.If it were possible to identify such a cut, then we would conclude that (18) does not contain an enumerated set of extreme points to λ-(11) for λ ∈ [0, 1] or that (θ 1 , θ 2 , y) is not an efficient extreme solution.As the BBSA is guaranteed to move from a current explored efficient solution to a neighbouring efficient solution of (18)-{ Dr , Λ, Dp } following the pivoting steps of the bi-objective parametric simplex method, we conclude that the BBSA cannot "miss" an efficient extreme solution as it proceeds and Lemma 1 confirms that every solution found is efficient.To see this we briefly outline two cases.In Case 1, a particular non-dominated extreme point ẑ (and corresponding efficient solution) is skipped, and its neighbouring points are correctly identified.This could be due to a cut that removes ẑ.However, this cut cannot exist (as also outlined in the second paragraph of the proof of Lemma 1).As we also know the parametric simplex method will correctly visit all non-dominated extreme points, ẑ cannot be missed.In Case 2, non-dominated point ẑ is missing because it is dominated by the objective vector z * of another solution that is incorrectly identified as efficient by the algorithm -this contradicts Lemma 1, and would indicate that cuts exist that would remove z * , hence z * could not have been confirmed as non-dominated.The cases are also illustrated in Figure C.17 in Appendix C.
The BBSA therefore finds a complete set of efficient extreme solutions with images corresponding to non-dominated extreme points of BLP (9).Proposition 3. The BBSA Algorithm 1 applied to BLP (9) terminates in a finite number of steps.
Proof.We start by making the following observations: 1.At each iteration of the BBSA the master problem BM defined by (18)-{ Dr , Λ, Dp } has a finite set of efficient extreme solutions, and so does each subproblem BS, if it is feasible.
2. Any feasible subproblem BS has a finite number of efficient extreme solutions from which a finite number of weighted optimality cuts is generated.
3. Both BM defined by (18)-{ Dr , Λ, Dp } and BS (11) have a finite set of efficient extreme solutions that can be obtained by the bi-objective parametric simplex method in a finite number of steps (or basis pivots) assuming appropriate strategies to prevent cycling are employed.
Initially, the BBSA solves a single-objective Benders problem λ-(18)-{∅, ∅, ∅} in a finite number of steps by adding a finite number of cuts.At each iteration the BBSA solve of BM (18)-{ Dr , Λ, Dp } in lines 23-31 starts from the most recently confirmed efficient solution ( θ1 , θ2 , ȳ) with image z, which is an optimal solution to λ-BM.The bi-objective simplex method performs one or more basis pivots thereby iterating, in a finite number of steps, to an unexplored efficient solution of BM (18)-{ Dr , Λ, Dp }, ( θℓ 1 , θℓ 2 , ȳℓ ) with image z ℓ , associated weight λ ′ ≦ λ, and strictly better second objective value (if it exists; otherwise the algorithm terminates).The corresponding efficient solution ( θℓ 1 , θℓ 2 , ȳℓ ) is then tested for feasibility and efficiency at the start of the next iteration of the BBSA, and two cases may occur after cuts are generated in lines 5-13.
Case 1: If solution ( θℓ 1 , θℓ 2 , ȳℓ ) remains feasible, and the algorithm successfully confirms this solution as efficient, then a new efficient solution with lower second objective value compared to ( θ1 , θ2 , ȳ) is obtained.
Case 2: If solution ( θℓ 1 , θℓ 2 , ȳℓ ) is neither feasible nor efficient, a violated feasibility cut or a finite number of weighted optimality cuts will be found to remove the solution from further consideration.Updating the master problem BM (18)-{ Dr , Λ, Dp } with these new cuts will yield a new candidate unexplored solution with weight λ ′′ ≦ λ ′ , while the last explored efficient solution ( θ1 , θ2 , ȳ) remains feasible.
In Case 2, the algorithm returns to this last explored efficient solution ( θ1 , θ2 , ȳ) of BM (18)-{ Dr , Λ, Dp } and performs a finite number of bi-objective parametric simplex iterations again to identify another candidate efficient solution in a finite number of steps (if it exists), as noted above.Case 2 will only occur a finite number of times until a neighbouring efficient extreme solution of ( θ1 , θ2 , ȳ) is reached.This is because, firstly, there exists only a finite number of feasibility cuts.Secondly, only a finite number of weighted optimality cuts are added for a particular weight λ ′ : The weight λ ′ corresponds to optimal solutions of the weighted sum problem for points on the facet between images z and zℓ of ( θ1 , θ2 , ȳ) and ( θℓ 1 , θℓ 2 , ȳℓ ) in objective space.This weight λ ′ may remain unchanged between multiple subsequent iterations of Case 2.
A finite set of weighted optimality cuts are added after solving BS that, in particular, include cuts for weight λ ′ .Therefore, single-objective Benders theory, see Benders (1962), confirms only a finite number of cuts will be added here by solving a finite number of subproblems to confirm λ ′ is the facet's true weight in the original problem, hereby confirming solution ( θ1 , θ2 , ȳ) is optimal for λ ∈ [λ ′ , λ], and that ( θℓ 1 , θℓ 2 , ȳℓ ) is efficient.Alternatively, a cut alters the facet and its associated weight to λ ′′ < λ ′ and the same argument can be made for this next λ ′′ .In summary, to identify each (of the finitely many) efficient basic feasible solution in Xe the BBSA conducts a finite number of steps in Case 2 to eventually arrive at Case 1 thus concluding efficiency of a solution in a finite number of iterations.

Computational Experiments
Benders decomposition has been successfully applied to a wide range of single-objective optimisation problems.We test the BBSA on three different types of BLP with decomposable structure.Sets of problem instances of varying size are tested for bi-objective variants of the Fixed Charge Transportation Problem (FCTP) (e.g.Zetina et al., 2019), the Portfolio Problem (PP) (e.g.Elci, 2022) and the Multiple Knapsack Problem (MKP) (e.g.Angulo et al., 2016;Maher, 2021).Extensive numerical testing is conducted to analyse the performance of the proposed BBSA.We are primarily concerned with solving BLPs at this stage and, as such, consider continuous relaxations of our problem instances in the testing.We compare the BBSA to a dichotomic approach that is referred to as Dichotomic Benders Algorithm (DBA) where a series of single-objective weighted-sum problems are solved with Benders decomposition.The numerical experiments were carried out in MATLAB R2020b and performed on a 64-bit Windows 10 PC system with Intel Xeon W-2145 CPU 3.7 GHz processor and 32 GB RAM.The open-source lpsolve package is used whenever we need to solve weighted BM and BS optimisation problems.
An implementation of the bi-objective simplex algorithm is used to find a set of non-dominated extreme points in BM and BS.
Before we discuss results in Section 4, the problem instances are briefly introduced in Sections 4.1-4.3.All instances are available in https://bitbucket.org/araith/bbsa-instances.The results of the BBSA for the different instance types are discussed in Section 4.4-4.6.Finally, the comparison of runtime between the BBSA and the DBA is presented in Section 4.8.

Bi-objective Fixed Charge Transportation Problem (Bi-FCTP)
In the FCTP it is assumed that there are m sources (e.g, factories) and n sinks (e.g., customers).
Each source can transfer commodities to any of the sinks with a per-unit transportation cost and a fixed charge.A solution to the FCTP involves satisfying the demand of each sink given the available capacity of each source.The first objective of the Bi-FCTP considered here minimises the total transportation cost (including variable and fixed costs) and the second one minimises the transportation time, where the fixed-cost component of the time objective is considered to be the preparation time required to set up the connection.The full formulation is shown in Appendix D.1.We adapt instances from the literature: FCTP Data sets 1 and 2 are from Agarwal and Aneja (2012), and data set 3 is from Roberti et al. (2015).Data set 1 contains 12 instances that have between 4 and 17 sources and between 3 and 43 sinks.Data sets 2 and 3 each have 30 instances, with n = m = 15 and n = m = 20, respectively.The unit costs are scaled to approximately maintain a pre-defined ratio σ between the total variable and fixed costs.It should be noted that this ratio was denoted θ in the original paper, but we use σ here to avoid confusion with our use of variables θ 1 , θ 2 in Section 2.3.Three different scaling factors are considered, σ = 0, 0.2, 0.5.Instances with σ = 0 are pure FCTP in which all unit costs are set to zero.For σ = 0.2 (σ = 0.5), the unit costs are scaled so that their total constitutes 20% (50%) of the total costs.The second objective function for Bi-FCTP is randomly generated with values in the same range as the first objective.

Bi-objective Portfolio Problem (Bi-PP)
In the PP, there are n potential assets, and the problem is how to distribute a unit investment among the n assets to minimise cost.The original PP model (Qiu et al., 2014)

Bi-objective Multiple Knapsack Problem (Bi-MKP)
The MKP is the third type of benchmark instance and is based on Angulo et al. (2016) who formulate a stochastic binary version of the MKP.In the problem, three sets of items are packed into two sets of knapsacks with the aim of minimising overall cost, and with separate weight constraints associated with each knapsack.Here we solve a relaxation of the MKP in the deterministic case.Angulo et al. (2016) generated 30 instances of identical problem dimensions, which are available in Ahmed et al. (2015).We follow the same approach as for Bi-FCTP and randomly generate coefficients for the second objective function of the Bi-MKP instances.However, in the original instances, the subproblem objective function coefficients are all one.Instead, we randomly generated the second objective function coefficients with uniform distribution over a larger range of values.The model is shown in Appendix D.3.

Bi-FCTP Results
An overview of average statistics within each instance group is given in Table 1.The data set is shown in column DS, the instance group name is listed in column n m.Instance groups are considered for instances 15 15 and 20 20 with ten instances for each value of σ, and their averages are reported in Table 1, whereas for the instances in data set 1 only a single instance exists.The table lists the number of non-dominated extreme points |Z n |, and the number of feasibility, optimality and active cuts (FC, OC and AC).Column "it" lists the number of iterations, tBBSA is the total runtime of the BBSA (in sec), and columns Ben%, BM% and BS% list the corresponding percentage of time taken by the initial single-objective Benders problem solves (Step 3 of Algorithm 1), BM and BS, respectively.The number of non-dominated extreme points found varies between 4 and 100 for the instances, where the average number of solutions in data sets 2 and 3 is fairly similar for different values of σ.Increasing problem size with larger number of sources n and sinks m leads to higher runtime, and so does an increase in σ as unit costs are only introduced into the problem when σ > 0 and increase as σ increases (Figure 4).We observe that there is no strong trend relating the number of non-dominated extreme points with the runtime, except in the instances where σ = 0.5.
Figure 5 shows the percentage of the runtime spent in different components of the BBSA algorithm for the different Bi-FCTP data sets.In data set 1 the results appear to vary; however, it should be noted that some of the instances are quite easy to solve, and thus run quickly (see Table 1).For the more challenging instances runtime in BM exceeds that in BS, often quite significantly.For data sets 2 and 3 and σ = 0, runtime of BM and the initial single-objective Benders solve are similar, and exceed the percentage of time spent in BS.The percentage of runtime spent in BM increases with increasing σ as the number of optimality and feasibility cuts tends to increase with increasing master problem complexity.
Only a small portion of optimality cuts are generated in the initial single-objective Benders problem, with the exception of the instances with σ = 0 in data sets 2 and 3. Here, two optimality cuts are generated in the initialisation and no further optimality cuts are generated as objective functions coefficients of BS are zero.For data sets 2 and 3, on average 73.88% and 78.07% of feasibility cuts are generated in the initialisation, whereas only 45.52%, on average, are generated for the instances in data set 1.For data set 3 with σ = 0 most feasibility cuts are generated in the initialisation (98.65% on average, whereas it is only 86.50% for data set 2).
In the BBSA iterations, for data set 1, feasibility cuts are added in 52.5% of all iterations, and 2.11 optimality cuts are generated per iteration in which optimality cuts are added.For data sets 2 and 3 feasibility cuts are added in 38.83% and 30.92% of all iterations.No optimality cuts are added for instances with σ = 0, whereas on average 2.33 and 2.71 optimality cuts are added per iteration for data sets 2 and 3 when σ > 0. The low number of optimality cuts per iteration is due to there being one or only a few non-dominated points in BS.A large number of feasibility cuts is required here since the total supply equals the total demand in the Bi-FCTP instances, making it difficult to find a feasible solution when decomposing the problem.For more challenging instances, such as data sets 2 and 3 with σ = 0.5 a large number of cuts is generated but only a fairly small number of active cuts is needed for the Benders reformulation of the bi-objective problems.There are between 31 and 378 active cuts in data sets 2 and 3 with an average of 2.85 cuts per non-dominated extreme point (5.05 for data set 1).
We can also compare the number of non-dominated points that were explored during the BBSA to the number of final non-dominated points, as identified at the end of Algorithm 1 (line 34), as shown in the ratio of explored over final non-dominated points in Figure 6.We see a clear increase in this ratio based on problem complexity in Bi-FCTP data sets 2 and 3 for higher values of σ indicating increasing problem complexity.
The considered problem instances suffer from degeneracy in BM and BS. Figure 7 shows the average number of degenerate pivots in BM per BBSA iteration, and the average number of degenerate pivots per BS solved.BS generally has few non-dominated extreme points.Due to the degeneracy, many iterations are, however, needed to solve BS with the bi-objective simplex algorithm.In a degenerate pivot, a variable enters the basis at value zero, and the objective values do not change.To confirm all non-dominated extreme points, hundreds of degenerate pivots may be required, and this is time-consuming.

Bi-PP Results
There are 100 instances in this class, 10 for each value of k.Table 2 reports averaged statistics for each instance group.Bi-PP k indicates the k-th group of Bi-PP instances where the total amount by which the return constraints can be violated is k.The number of non-dominated extreme points varies from 32 to 90, with an average of 59.98, which shows a fair spread in the number of non-dominated extreme points.As k increases, the number of non-dominated extreme points also tends to grow.As expected, with increasing number of non-dominated extreme points, the time taken to solve the problems also rises, as demonstrated in Figure 8.
However, the difference between the average time taken to solve instances with k = 1 and instances k = 10 is not large.
The percentage of time taken in different components of the BBSA for Bi-PP is depicted in Figure  problem is more time consuming when k is small, as shown in Figure 9.A.
Figure 9.B shows the average number of explored non-dominated points per extreme point for different groups of instances.We observe that the number of explored non-dominated points per non-dominated extreme point is spread from 3.17 to 14.93 with a mean of 12.36.The average number of non-dominated points that are explored to identify an non-dominated extreme point increases slightly with increasing k.
It should also be noted that there is a single optimal solution in BS for all Bi-PP instances.
When the values of the y variables in BM are fixed, the values of the x variables can be uniquely derived and a single non-dominated extreme point is identified in BS.Therefore, two optimality cuts are generated by BS in each iteration corresponding to θ 1 and θ 2 with λ = 0 and λ = 1.
Bi-PP instances are therefore examples of Remark 2, where the optimal solution corresponding  to λ = 1 and λ = 0 of BS are the same (i.e.there is only one non-dominated point in BS).
Therefore, to identify a set of non-dominated extreme points, it is sufficient to generate single optimality cuts for θ 1 and θ 2 from BS in each iteration.
Figure 10.A gives an overview of the number of cuts generated for Bi-PP instances.An interesting observation is that when k is small, more feasibility cuts are generated as it is more likely to encounter infeasibility when the permitted constraint violation in BS is small.However, when k is larger, the number of feasibility cuts generated from BS drops significantly.There is a large variation of between 1.78 and 22.34 active cuts, with an average of 7.99 per non-dominated extreme point.More initial optimality and feasibility cuts are generated in Bi-PP instances in comparison with Bi-FCTP.18.09% of the total cuts are generated in the initial single-objective Benders solve (Algorithm 1, line 2), where 99.04% and 13.47% of the feasibility and optimality cuts are generated initially.
In contrast with Bi-FCTP, Bi-PP instances have fewer degenerate pivots as shown in Figure 10.B for BM where the number of degenerate pivots per BBSA iteration is between 2.39 and 4.79, with a mean of 3.41, which is significantly less than for Bi-FCTP instances.The BS has one non-dominated point, and there is no pivoting required to confirm this, and as such, there are no degenerate pivots in BS.

Bi-MKP Results
Table 3 summarises the numerical results of Bi-MKP.We report each column's minimum, maximum and mean value.For Bi-MKP instances, the number of non-dominated extreme points (column |Z n | in Table 3) is much higher than for the Bi-FCTP and the Bi-PP instances.There  Figure 11.A gives an overview of cut generation for Bi-MKP instances.There are between 23.11% and 39.89% active cuts, with an average of 30.81%.No feasibility cuts need to be generated for these instances.There are a large number of optimality cuts, only about 0.006% of which are generated when solving the initial single-objective Benders problem.
The number of explored non-dominated points per non-dominated extreme point is much lower compared to the other instance types.For Bi-MKP instances, only between 1.34 and 1.63 nondominated points are required to be explored to identify an non-dominated extreme point, with an average of 1.47.
Compared to the Bi-PP instances, there are even fewer degenerate pivots in these instances.
From Figures 11.B and 11.C, it is clear that BS has a similar number of degenerate pivots per BBSA iteration as the BM.The number of degenerate pivots per BBSA iteration in BS is between 0.3 and 1.2, with an average of 0.58.For BM, the number of degenerate pivots per iteration is between 0.29 and 1.36, with an average of 0.57.These results show that, when the number of degenerate pivots in BS and BM is small enough, the BBSA quickly identifies a set of non-dominated extreme points, even if there are many non-dominated extreme points.
Another interesting feature of the Bi-MKP instances is the number of weighted optimality cuts generated from BS.In contrast to the Bi-FCTP and Bi-PP in stances, a significant number of weighted optimality cuts are generated in Bi-MKP instances.For Bi-MKP instances, between 65.86 and 95.64, with an average of 78.68 weighted optimality cuts, are generated per iteration.

Summary of experiments for Bi-FCTP, Bi-PP and Bi-MKP instances
Several observations can be made from results obtained from the three different types of instances.For Bi-FCTP instances we observe that increasing problem size makes the problems more challenging to solve, and so does an increase of σ.Not surprisingly, the number of generated cuts generally increases for all problem types as the number of non-dominated points increases, although more feasibility cuts and active cuts are required for Bi-PP instances with lower k that tend to have fewer non-dominated points.It is interesting that BBSA runtime spent in BM and BS varies by instance type.Runtime spent in BM is higher for Bi-FCTP instances, and it is higher for BS in Bi-PP instances, whereas both are similar for Bi-MKP instances.Runtime performance appears to depend on how many points need to be explored by the BBSA compared to the number of final non-dominated extreme points.Runtime also depends on problem degeneracy where we observe that a high number of degenerate pivots are required for Bi-FCTP instances, whereas the other two instance types have much fewer degenerate pivots.In our experiments the highest runtimes are observed for the most challenging Bi-FCTP instances.BS has few efficient solutions for Bi-FCTP instances, often only one, meaning that few weighted optimality cuts are added in most iterations.In contrast, Bi-MKP instances have many efficient solutions in BS and therefore many weighted optimality cuts.

Comparison of the BBSA and DBA
We now compare the performance of the proposed BBSA with DBA.The idea of DBA is to iteratively solve weighted-sum problems (9) with the standard (single-objective) Benders decomposition technique where the weights of λ are identified from the dichotomic technique for BLPs (e.g.Aneja and Nair, 1979;Cohon et al., 1979) to iteratively identify all non-dominated extreme points.In the following, we compare the runtime (in seconds) of the proposed BBSA with the DBA on the Bi-FCTP, Bi-PP, and Bi-MKP instances.
Tables 4-6 give a summary of the DBA runs for our three instance types.Table 4 has the same first six columns as Table 1, and Tables 5 and 6 have the same first four columns as Tables 2   and 3.The subsequent columns are the number of DBA iterations (itDBA), that is the number of times single-objective Benders decomposition was called to solve a weighted problem (note that itDBA = 2|Z n | − 1); the DBA runtime (tDBA); the ratio of DBA and BBSA runtime (tDBA/tBBSA) and the average DBA runtime per iteration (tDBA per it).If the runtime ratio tDBA/tBBSA exceeds 1, then BBSA runs faster than DBA.Tables 5 and 6 show averages.
Solving our problem instances with the proposed BBSA is more effective than iteratively solving single-objective weighted-sum standard Benders decomposition problems, as the columns tDBA/tBBSA in Tables 4-6 show very clearly, where the run-time ratio mostly exceeds 1, often quite significantly.Some Bi-FCTP instances are challenging and a time limit of 10000 seconds is applied to DBA.The DBA did not find all non-dominated extreme points in two instances of Bi-FCTP data set 3. There are only a few small instances where DBA outperforms the BBSA.
Depending on the type of instance, different factors affect DBA performance.For Bi-FCTP instances, degeneracy makes solving problems with Benders decomposition slow in general, and this applies also to the single-objective weighted Benders problems solved by DBA.We also  observe that DBA generates many more feasibility and optimality cuts than the BBSA.In making this comparison it should be noted that different iterations of DBA may generate the same cut and we might hence be counting some cuts multiple times, whereas cuts counted for the BBSA are always unique cuts.For instance, in Data Set 2 of the Bi-FCTP instances, between 21.5 and 42.26 times as many feasibility cuts are generated by DBA compared to the BBSA, and this ratio is significantly higher for Data Set 3. DBA also generates more optimality cuts although the ratio is lower (between 2.25 and 3.93 for instances with σ > 0).Both approaches would be similarly affected by degeneracy.
For Bi-PP to find a non-dominated extreme point, lots of standard Benders decomposition iterations are required, particularly when k is small as can be seen in Figure 9.A and Table 2 where the initial single-objective Benders problem solve makes up a large component of the  5.
DBA experiences large runtimes per iteration, especially in the more challenging instances with smaller values of k where large numbers of cuts are generated across the iterations (Table 5).
Finally, for the Bi-MKP instances, we observe that DBA generates significantly fewer optimality cuts when compared with BBSA -only between 10.17% and 16.13%.DBA also generates under 10 optimality cuts per iteration, hence converging to a solution for each single-objective problem quickly.Despite generating fewer cuts overall (and in each iteration), DBA does show longer overall runtime.This is because the problem has many non-dominated extreme points.If there are p non-dominated extreme points, DBA needs to solve almost twice as many (2p−1) weighted single-objective problems with Benders decomposition.Even though individual problems solve fast, runtime accumulates here.Furthermore, the BBSA is implemented to only retain currently active cuts to keep the size of the master problem BM small (hence enabling us to solve it quickly in each iteration), and the subproblems BS also solve very quickly due to low degeneracy.

Conclusion and Future Work
In this paper, we have developed the necessary theory for applying Benders decomposition in a bi-objective setting.In particular, we present the Benders reformulation of a generic BLP and argue that a finite subset of the constraints is sufficient to generate a complete set of efficient extreme solutions.Furthermore, we discuss the need for so-called weighted optimality cuts and propose an iterative procedure that combines the bi-objective simplex algorithm with the cut generation procedure of the standard Benders decomposition algorithm.We have shown that the proposed methodology can obtain a complete set of efficient extreme solutions and the corresponding set of non-dominated extreme points to a BLP.Using sets of diverse problem instances, we have also identified factors that affect the performance of the proposed algorithm and showed that it performs well when compared to a dichotomic weighting approach combined with standard (single-objective) Benders decomposition.
There are several topics that can be explored in future research.Firstly, the results of BBSA have shown that, as problem size increases, runtime of the BBSA can increase quite quickly.
It could therefore be beneficial to investigate ways to improve the performance of the Benders decomposition approach.Acceleration techniques for the single-objective case can be adapted to the bi-objective setting.Secondly, methods that exploit the bi-objective nature of the problem are also of interest.A bi-directional approach that simultaneously explores non-dominated points from left to right and right to left is one such possibility.So too are methods that can improve the accuracy of λ weights considered.The bi-objective simplex algorithm determines the next λ weight for the problem to which it is being applied.Investigating whether or not one can accurately determine the next λ weight within a Benders decomposition is certainly interesting and would remove some of the redundant steps in the BBSA.Knowing the next λ would enable one to simply iteratively solve single-objective weighted-sum problems to obtain a complete set of efficient extreme solutions.From an application perspective, it could be interesting to apply the developed methodology to stochastic BLPs with multiple subproblems.It would also be of interest to extend the proposed approach to multi-objective LPs with three or more objectives.As long as an algorithm to solve such multi-objective LPs is available, such as the multi-objective simplex algorithm (Evans and Steuer, 1973;Rudloff et al., 2017) one should be able to develop a multi-objective version of Algorithm 1 to iteratively identify efficient solutions in a multi-objective master problem, generate cuts in a multi-objective subproblem and, by doing so, iteratively identify a complete set of efficient solutions.These multi-objective simplex algorithms are more challenging to implement because efficient solutions of a multi-objective optimisation LP have multiple neighbouring efficient solutions that all need to be explored.

Figure 2 :
Figure 2: Feasible region of dual formulation for the first subproblem (left) and the second (right).
) therefore includes Benders feasibility cuts of the form (b − B ȳ) ⊤ π r ≦ 0 as well as optimality cuts of the form −θ + (b − B ȳ) ⊤ π λ p ≦ 0 for extreme rays π r ∈ D r and extreme points π λ p ∈ D λ p of the dual of the λ-weighted subproblem.
considers an uncertain return, where the overall return should be at least r.Here we consider the case where all m possible return scenarios are included simultaneously.The formulation allows violation of some of the overall return constraints, and the total amount of violation is limited to be at most k.The full model is shown in Appendix D.2.Since the original problem has a single-objective function, we randomly create a second one to obtain Bi-PPs.In Qiu et al. (2014), 10 different instances are generated where m = 200, and n = 20 and the number of the violated constraints is k = 15.We generated 100 instances, where the total amount by which the return constraints are violated ranges from k = 1 to k = 10 with 10 different instances per k.

Figure 6 :
Figure 6: Ratio of explored non-dominated (ND) over final non-dominated extreme points.

Figure 7 :
Figure 7: The average number of degenerate pivots per iteration of the BBSA in which BM is solved is shown in the top row; The average number of degenerate pivots in BS in the bottom row.

Figure 8 :
Figure 8: Number of non-dominated extreme points and runtime of Bi-PP instances.

Figure 9 :
Figure 9: Comparison of the BBSA components runtime for Bi-PP instances (A).The number of explored nondominated points per extreme points for Bi-PP instances (B).
is no clear trend in the relationship between runtime and |Z n |.The percentage of the runtime spent in different components of the BBSA is different from other problem instance types.For Bi-MKP instances, solving BM and BS takes a similar amount of time, and the time to solve the initial single-objective Benders problems is very low.

Figure 10 :
Figure 10: The number of cuts of the BBSA for Bi-PP instances (A).The number of degenerate pivots of BM for Bi-PP instances (B).

Figure A. 12
Figure A.12: A geometric illustration

Table 1 :
Summary of run statistics for Bi-FCTP instance groups.Averages are shown for instance groups 15 15 and 20 20 over sets of 10 instances for each value of σ.

Table 2 :
Summary of run statistics for Bi-PP instance groups, where the average values are shown.

Table 3 :
Summary of minimum, mean and maximum run statistics for Bi-MKP instances.

Table 4 :
Summary of DBA run statistics for Bi-FCTP instance groups.Averages are shown for instance groups 15 15 and 20 20 over sets of 10 instances for each value of σ.

Table 5 :
Summary of DBA run statistics for Bi-PP instance groups, where the average values are shown.

Table 6 :
Summary of minimum, mean and maximum DBA run statistics and Bi-MKP instances.for small k, as well as in the time of DBA per iteration (tDBA per it) in Table