Given λ, S L , and S U , the interval representation of f ( x P opt ( λ ) ) is R ( S L , S U , λ ) = ( [ L 1 ( S L , λ ) , U 1 ( S U , λ ) ] , … , [ L k ( S L , λ ) , U k ( S U , λ ) ] ).

For k = 2, components of the interval representation of f ( x P opt ( λ ) ), lower and upper shells, as well as vectors of lower and upper bounds, are illustrated in Fig. 1.

To gauge the quality of R ( S L , S U , λ ), we calculate (5) G P s u b , l ( λ ) : = 100 × U l ( S U , λ ) − L l ( S L , λ ) U l ( S U , λ ) , l = 1 , … , k .

G P s u b ( R ( S L , S U , λ ) ) : = ( G P s u b , 1 ( λ ) , … , G P s u b , k ( λ ) ) forms the Pareto suboptimality gap of interval representation R ( S L , S U , λ ).

As in Kaliszewski and Miroforidis (2022), we will use S L : = { I N C λ } as a valid lower shell one can use to calculate L l ( S L , λ ), l = 1 , … , k.

To populate upper shell S U , the following two lemmas defined in Kaliszewski and Miroforidis (2022) can be used.

Given X 0 ′ ⊃ X 0 and λ, let ChebRLX ( X 0 ′ , λ ) denote the Chebyshev scalarization (problem (2)) of some relaxation of the MOMIP problem (problem (1)) with feasible set X 0 ′ for some λ. Based on Lemmas 2 and 3, one can derive a single-element upper shell by solving ChebRLX ( X 0 ′ , λ ). Given X 0 ′ ⊃ X 0 , the sum of such single-element upper shells derived for different vectors λ forms upper shell S U .

The surrogate relaxation of problem (1) with X 0 ′ ( μ ) : = { x ∈ X | ∑ p = 1 m μ p g p ( x ) ⩽ ∑ p = 1 m μ p b p } instead of X 0 , where μ p ⩾ 0, p = 1 , … , m, μ ≠ 0, is a valid relaxation of this problem (see Kaliszewski and Miroforidis, 2022). μ is a vector of surrogate multipliers. We will use this type of relaxation with μ as a parameter to derive elements of an upper shell. We also follow what has been shown in Kaliszewski and Miroforidis (2022) on large-scale instances of the MOMIP problem that for a given μ solving the Chebyshev scalarization of the surrogate relaxation of the MOMIP problem by a MIP solver is much easier than solving the Chebyshev scalarization of the MOMIP problem.

Given λ and S L , based on Lemma 1, element x of upper shell S U is a source of an upper bound on f l ¯ ( x P opt ( λ ) ), l ¯ ∈ { 1 , … , k }, when f ( x ) is appropriately located with respect to the vector of lower bounds L ( S L , λ ). In Miroforidis (2021), an idea of how to derive an upper shell that consists of an element useful to calculate upper bounds on f l ¯ ( x P opt ( λ ) ) has been proposed. This idea is to probe the objective space by perturbing components of vector λ. Yet, there is no algorithmic approach in Miroforidis (2021) doing that. In the current work, we try to fill in this gap.

For k = 2, the idea of an algorithm for deriving upper shell S U 1 whose some element is a source of an upper bound on f 1 ( x P opt ( λ ) ) is shown in Fig. 2. Let us assume that vector μ is given. At the beginning, S U 1 : = ∅. Below, we describe three iterations of this example.

Iteration 1. We set the first probing vector λ ′ : = ( λ 1 + δ , λ 2 − δ ), λ 2 ′ > 0, for some δ > 0 (as a probing vector, we exclude λ because we expect that the corresponding solution will not be properly located with respect to the vector of lower bounds, see Kaliszewski and Miroforidis, 2021). Let x λ ′ be the solution to ChebRLX ( X 0 ′ ( μ ) , λ ′ ). S U 1 : = S U 1 ∪ { x λ ′ }. Based on Lemma 1, x λ ′ is not a source of an upper bound on f 1 ( x P opt ( λ ) ) because f ( x λ ′ ) is not appropriately located with respect to the vector of lower bounds L = L ( S L , λ ). Hence, we CONTINUE.

Iteration 2. We set λ ″ : = ( λ 1 ′ + δ , λ 2 ′ − δ ), λ 2 ″ > 0. Let x λ ″ be the solution to ChebRLX ( X 0 ′ ( μ ) , λ ″ ). S U 1 : = S U 1 ∪ { x λ ″ }. Based on Lemma 1, x λ ″ is not a source of an upper bound on f 1 ( x P opt ( λ ) ). Hence, we CONTINUE.

Iteration 3. We set λ ‴ : = ( λ 1 ″ + δ , λ 2 ″ − δ ), λ 2 ‴ > 0. Let x λ ‴ be the solution to ChebRLX ( X 0 ′ ( μ ) , λ ‴ ). S U 1 : = S U 1 ∪ { x λ ‴ }. Based on Lemma 1, x λ ‴ is a source of an upper bound on f 1 ( x P opt ( λ ) ). So, U 1 = U ( S U 1 , λ ) = f 1 ( x λ ‴ ). Hence, we STOP.

As elements x λ ′ , x λ ″ , and x λ ‴ are Pareto optimal solutions to the relaxation of the MOMIP problem with X 0 ′ ( μ ), S U 1 is a valid upper shell.

To obtain an upper bound on f 2 ( x P opt ( λ ) ), we need to derive upper shell S U 2 . To do this, in the first iteration, we set λ ′ : = ( λ 1 − δ , λ 2 + δ ), λ 1 ′ > 0, and proceed in the same way.

Given l ¯ ∈ { 1 , … , k }, the FindUpperShell algorithm tries to derive upper shell S U l ¯ whose element x is a source of an upper bound on f l ¯ ( x P opt ( λ ) ), i.e. U l ¯ ( S U l ¯ , λ ).

In Line 2, we set step size δ that is used to modify components of consecutive probing vectors λ ′ . Parameter γ > 0 controls the step size, i.e. the greater the value of parameter γ, the denser the sampling of the objective space to search for the desired element of the upper shell. In the main loop (Lines 4–16), we populate upper shell S U l ¯ checking if its new element fulfills conditions of Lemma 1 to be a valid source for the upper bound on f l ¯ ( x P opt ( λ ) ). The algorithm stops when λ l ¯ ′ ⩾ 1 OR some component of λ ′ is negative OR an element that fulfills conditions of Lemma 1 is found. Lines 6 and 9 guarantee that ∑ l = 1 k λ l ′ = 1. The exit condition of the “while” loop ensures that λ l ′ > 0, l = 1 , … , k.

If no element of S U l ¯ satisfies conditions of Lemma 1, then the algorithm simply returns the upper shell as well as the only available upper bound on f l ¯ ( x P opt ( λ ) ), namely y l ¯ ∗ . Otherwise, an upper bound better than y l ¯ ∗ is returned, as well as the upper shell.

Based on the above elements, an interval representation of the Pareto optimal outcome given by vector λ can be calculated with the use of the Chute algorithm. Along with the interval representation, this algorithm also returns lower and upper shells that were determined during its operation. In Section 5.7, we explain why the algorithm also returns lower and upper shells.

Line 4 of the algorithm needs clarification. Vector μ can be set as shown in Kaliszewski and Miroforidis (2022), namely by taking μ p : = 1, p = 1 , … , m. In that work, all surrogate multipliers have the same value. We call this version of the Chute algorithm Chute1.

Yet, in Kaliszewski and Miroforidis (2022), in the section “Final remarks”, it has been suggested that “Tighter bounds might be obtained with other values of the multipliers. This possibility is worth exploring in future works.”. Unfortunately, there is no idea there how to select a vector of surrogate multipliers other than ( 1 , … , 1 ) ∈ R m . However, we can use the theory of duality for this purpose.

Given μ, λ, let x be the solution to ChebRLX ( X 0 ′ ( μ ) , λ ), and s ( μ ) be the objective function value of x. Based on Lemmas 2 and 3, { x } is a valid upper shell. Let s be the objective function value of the solution to problem (2) with λ and X 0 . It is a well-known fact (see, e.g. Glover, 1965, 1968) that s ⩾ s ( μ ). Hence, for a given λ and μ, s ( μ ) is a lower bound on values of s.

Given λ, the best (highest) lower bound s ∗ on values of s is the objective function value of the solution μ ∗ to the following surrogate dual problem (6) sup μ ⩾ 0 , μ ≠ 0 { min x ∈ X 0 ′ ( μ ) max l λ l ( y l ∗ − f l ( x ) ) + ρ e k ( y ∗ − f ( x ) ) } that is connected to the Chebyshev scalarization (problem (2)). Solving (6) to optimality can be time-consuming. Yet, a suboptimal vector of multipliers μ ˜ can be determined instead of μ ∗ . It can be done with the help of a quasi-subgradient-like algorithm (we shall call it Suboptimal) by Dyer (1980) with the following stopping condition.

“Number of iterations without improving the value of the objective function in problem (6) is greater than N” OR “time limit on optimization is greater than T S seconds”.

In the current work, we set time limits on computation, hence the above stopping condition is justified in practice.

We will use vector ( 1 , … , 1 ) | | ( 1 , … , 1 ) | | ∈ R m , where | | . | | is the Euclidean norm, as an initial vector of surrogate multipliers in the Suboptimal algorithm. Under the above assumptions, this algorithm has three parameters, namely λ, N, and T S , i.e. in pseudocode, it can be used as a function Suboptimal ( λ , N , T S ), which returns a suboptimal vector of multipliers μ ˜.

We shall call a version of the Chute algorithm that uses (in Line 4) the Suboptimal algorithm to set vector of surrogate multipiers μ for a given λ Chute2. It has two additional input parameters N and T S .

Let us note that in the Chute2 algorithm, we set the vector of surrogate multipliers once for a given λ. The FindUpperShell algorithm uses perturbations of the λ vector to sample the objective space, and for all these perturbations the same vector μ is used. It is our heuristic assumption that even using the same vector μ for various vectors λ ′ , that are close to λ, the Chute2 algorithm is able to find a better R ( S L , S U , λ ), so tighter upper bounds on components of f ( x P opt ( λ ) ), than the Chute1 algorithm. However, it is at the cost of increasing the computation time relative to Chute1 by at most T S . We will check it experimentally in the next section.

The idea (for k = 2 and ρ = 0) of using suboptimal values of surrogate multipliers to get an upper shell that is a source of a better upper bound is illustrated in Fig. 3. Let μ I : = ( 1 , … , 1 ) ∈ R m , and μ ˜ be the output of the Suboptimal algorithm (with some N and T S ) for some λ. x ( μ I ) is the solution to optimization problem (2) with X 0 replaced with X 0 ′ ( μ I ), and s ( μ I ) is the objective function value of x ( μ I ). x ( μ ˜ ) is the solution to optimization problem (2) with X 0 replaced with X 0 ′ ( μ ˜ ), and s ( μ ˜ ) is the objective function value of x ( μ ˜ ). Vector of lower bounds L is marked with a triangle. Both elements f ( x ( μ ˜ ) ) and f ( x ( μ I ) ) are appropriately located with regards to L (see Lemma 1) to be sources for an upper bound on f 2 ( x P opt ( λ ) ). As s ( μ ˜ ) > s ( μ I ) (contours of the Chebyshev metric for both values s ( μ ˜ ) and s ( μ I ) are shown by solid thin lines) as well as f 1 ( x ( μ ˜ ) ) < f 1 ( x ( μ I ) ) AND f 2 ( x ( μ ˜ ) ) < f 2 ( x ( μ I ) ), element x ( μ ˜ ) is a source of a better upper bound on f 2 ( x P opt ( λ ) ) than element x ( μ I ), as upper bounds are calculated with the use of components of the upper shell elements. In Fig. 3, f ( x ( μ ˜ ) ) is closer to the (unknown) Pareto front (represented by the solid curve) than element f ( x ( μ I ) ). It could happen that condition f 1 ( x ( μ ˜ ) ) < f 1 ( x ( μ I ) ) AND f 2 ( x ( μ ˜ ) ) < f 2 ( x ( μ I ) ) does not hold. In this case, if f 2 ( x ( μ ˜ ) ) ⩾ f 2 ( x ( μ I ) ), we obtain no better upper bound on f 2 ( x P opt ( λ ) ). On the other hand, if f 1 ( x ( μ ˜ ) ) ⩾ f 1 ( x ( μ I ) ) and still f ( x ( μ ˜ ) ) is appropriately located with regards to L (see Lemma 1) to be a source for an upper bound on f 2 ( x P opt ( λ ) ), we get a better upper bound on f 2 ( x P opt ( λ ) ).

For k > 1, the MOMKP is formulated in the following way. (7) vmax f 1 ( x ) : = ∑ j = 1 n c 1 , j x j … f k ( x ) : = ∑ j = 1 n c k , j x j s.t. x ∈ X 0 : = { x | ∑ j = 1 n a p , j x j ⩽ b p , p = 1 , … , m , x j ∈ { 0 , 1 } , j = 1 , … , n } , where all a p , j , c p , j are non-negative. In Kaliszewski and Miroforidis (2022), it has been explained why the MOMKP is N P-hard.

As tri-criteria instances of the MOMKP, we take two instances from Kaliszewski and Miroforidis (2022) that were generated based on the 1st problem of the 6th group ( n = 500, m = 10) of multidimensional 0–1 knapsack problems, as well as on the 1st problem of the 9th group ( n = 500, m = 30) (both single-objective problems are stored in Beasley OR-Library, http://people.brunel.ac.uk/~mastjjb/jeb/info.html). We call these tri-criteria instances Three6.1 and Three9.1, respectively.

By removing the third objective function of problem Three6.1, we create a bi-criteria instance called Bi6.1. Analogously, by removing the third objective function of problem Three9.1, we create a bi-criteria instance called Bi9.1.

Bi6.1, Bi9.1, Three6.1, and Three9.1 are our test instances of the MOMIP problem.1¹

The instances can be made available to the reader by e-mail upon request.

Gurobi (version 10.0.0) for Microsoft Windows (x64) is our selected MIP solver. The optimizer is installed on the Intel Core i7-7700HQ-based laptop with 16 GB RAM.

To be consistent with limiting optimization time, we do not derive element y ˆ to optimality. So it is not known. Instead, we separately maximize each objective function of problem (7) within the time limit equal to 400 seconds. For instances Three6.1 and Three9.1, we set y l ∗ , l = 1 , 2 , 3, to the best upper bound (provided by the MIP solver) on values of the respective objective functions (none of these maxima the MIP solver determined in this time limit). Thus, for Three6.1, y ∗ : = ( 128872 , 131116 , 131738 ), and for Three9.1, y ∗ : = ( 119379.88 , 119365 , 118122 ). Obviously, y ˆ l < y l ∗ . Hence, such y ∗ approximating y ˆ can be used in (2) as well as when calculating lower and upper bounds. To avoid redundant calculations, for instances Bi6.1 and Bi9.1, we take only the first two components of respective y ∗ .

We set T L : = 1200 seconds, ρ : = 0.001. The absolute lower bound on values of all objective functions of the MOMKP is 0 (see Section 5.1), so this value is used for calculating lower bounds L l ( S L , λ ) , l = 1 , … , k (for details, see Kaliszewski and Miroforidis, 2022).

For instances Bi6.1 and Bi9.1, we generate one set of five vectors λ uniformly sampled from two-dimensional unit simplex (see Smith and Tromble, 2004), and obtain corresponding vectors of lower bounds L ( S L , λ ) (Table 1). For instances Three6.1 and Three9.1, we generate separate sets of five vectors λ, uniformly sampled from the three-dimensional unit simplex, and obtain corresponding vectors of lower bounds L ( S L , λ ) (Tables 2 and 3, respectively).

For all problem instances, for no vector λ, the selected MIP solver derived the solution to problem (2) in the assumed T L = 1200 seconds. Thus, the use of the Chute algorithm to determine interval representations of Pareto optimal outcomes given by vectors λ is justified.

We conduct two numerical experiments. In experiment 1, we test the behaviour of algorithm Chute1. In experiment 2, we test the behaviour of algorithm Chute2. In both experiments, on generated test instances, we test the behaviour of the Chute algorithm for γ : = 10 , 30 , 50. Given λ, we check the impact of parameter γ on components of G P s u b ( R ( S L , S U , λ ) ).

In the tables with results of the experiments, the meaning of the columns is as follows.

In bold, we indicate the improvement of a single component of G P s u b ( R ( S L , S U , λ ) ) when changing the value of parameter γ to a higher value, namely, we consider changes from γ = 10 to γ = 30 and from γ = 30 to γ = 50. By underscore, we indicate the deterioration of a single component of G P s u b ( R ( S L , S U , λ ) ) when changing the value of parameter γ to a higher value, namely, we consider changes from γ = 10 to γ = 30 and from γ = 30 to γ = 50.

In this experiment, we check the behaviour of the Chute1 algorithm.

The results for instance Bi6.1 are shown in Tables 4–6, and the results for instance Bi9.1 are shown in Tables 7–9. The results for instance Three6.1 are shown in Tables 10–12, and the results for instance Three9.1 are shown in Tables 13–15.

For instance Bi6.1, when changing γ = 10 to γ = 30, we observe an improvement in at least one component of G P s u b ( R ( S L , S U , λ ) ) for four vectors λ, although only in two cases (λ Nos. 3 and 5) two components improve. Yet, when changing γ = 30 to γ = 50, we observe an improvement of G P s u b ( R ( S L , S U , λ ) ) in at least one of its components for three vectors λ. For λ No. 3, we observe a deterioration of the first component of G P s u b ( R ( S L , S U , λ ) ), and at the same time, an improvement of the second. For λ No. 5, we observe an improvement of the first component of G P s u b ( R ( S L , S U , λ ) ), and at the same time, a deterioration of the second.

When changing γ = 10 to γ = 50, we observe an improvement in at least one component of G P s u b ( R ( S L , S U , λ ) ) for all vectors λ.

For instance Bi9.1, we observe a similar phenomenon (an improvement, as well as a deterioration), although when changing γ = 10 to γ = 30, we observe an improvement in only one component of G P s u b ( R ( S L , S U , λ ) ) for just two vectors λ. When changing γ = 10 to γ = 50, we observe an improvement in G P s u b ( R ( S L , S U , λ ) ) for vectors λ Nos. 1–4 (at least one component improves). For λ No. 5, G P s u b ( R ( S L , S U , λ ) ) remains unchanged.

For instances Bi6.1 and Bi9.1, the higher the value of parameter γ (higher sampling density of the objective space), the more numerous the derived upper shells are. For each vector λ, the time to derive the corresponding upper shell is a small fraction of the assumed time limit T L = 1200 seconds.

Let us check the results for tri-criteria instances. For instance Three6.1, when changing γ = 10 to γ = 30, we observe an improvement of G P s u b ( R ( S L , S U , λ ) ) in at least one of its components for three vectors λ, although only for one (λ No. 4) all components improve. When changing γ = 30 to γ = 50, we observe an improvement of G P s u b ( R ( S L , S U , λ ) ) in at least one of its components for three vectors λ. We also observe a deterioration of the first component of G P s u b ( R ( S L , S U , λ ) ) for λ No. 2, and, at the same time, an improvement on its third component. When changing γ = 10 to γ = 50, we observe an improvement in G P s u b ( R ( S L , S U , λ ) ) for all vectors λ (at least one component improves). For instance Three9.1 and vectors λ Nos. 2–5, upper bounds on all components of f ( x P opt ( λ ) ) are equal to the corresponding component of y ∗ . Only for λ No. 1, f 3 ( x P opt ( λ ) ) < y 3 ∗ , and when changing γ = 10 to γ = 30, as well as γ = 30 to γ = 50, there is an improvement only for the third component of G P s u b ( R ( S L , S U , λ ) ). For instances Three6.1 and Three9.1, the higher the value of parameter γ (higher sampling density of the objective space), the more numerous the derived upper shells are. For instance Three6.1 and instance Three9.1 with γ = 10, for most vectors λ, the time to derive the corresponding upper shell is a small fraction of the assumed time limit T L = 1200 seconds. However, for instance Three9.1 with γ = 30 and γ = 50, the time to derive the corresponding upper shell increases significantly compared to γ = 10, for all vectors λ.