A General Framework for Providing Interval Representations of Pareto Optimal Outcomes for Large-Scale Bi- and Tri-Criteria MIP Problems

Filcek, Grzegorz; Miroforidis, Janusz

doi:10.15388/24-INFOR549

A General Framework for Providing Interval Representations of Pareto Optimal Outcomes for Large-Scale Bi- and Tri-Criteria MIP Problems

Article type: Research Article

Authors: Filcek, Grzegorz¹ | Miroforidis, Janusz^{2; ∗}

Affiliations: [1] Wroclaw University of Science and Technology, Ul. Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland | [2] Systems Research Institute, Polish Academy of Sciences, Ul. Newelska 6, 01-447 Warszawa, Poland

Correspondence: [∗] Corresponding author.

Keywords: multi-objective mixed-integer programming, large-scale optimization, Chebyshev scalarization, Pareto front approximations, lower bounds, upper bounds

DOI: 10.15388/24-INFOR549

Journal: Informatica, vol. 35, no. 2, pp. 255-282, 2024

Received December 2023

Accepted March 2024

Published: 2024

Get PDF

Abstract

The Multi-Objective Mixed-Integer Programming (MOMIP) problem is one of the most challenging. To derive its Pareto optimal solutions one can use the well-known Chebyshev scalarization and Mixed-Integer Programming (MIP) solvers. However, for a large-scale instance of the MOMIP problem, its scalarization may not be solved to optimality, even by state-of-the-art optimization packages, within the time limit imposed on optimization. If a MIP solver cannot derive the optimal solution within the assumed time limit, it provides the optimality gap, which gauges the quality of the approximate solution. However, for the MOMIP case, no information is provided on the lower and upper bounds of the components of the Pareto optimal outcome. For the MOMIP problem with two and three objective functions, an algorithm is proposed to provide the so-called interval representation of the Pareto optimal outcome designated by the weighting vector when there is a time limit on solving the Chebyshev scalarization. Such interval representations can be used to navigate on the Pareto front. The results of several numerical experiments on selected large-scale instances of the multi-objective multidimensional 0–1 knapsack problem illustrate the proposed approach. The limitations and possible enhancements of the proposed method are also discussed.

1Introduction

The derivation of optimal solutions to large-scale instances of the Mixed-Integer Programming (MIP) problem can be impossible within a reasonable time limit even for contemporary commercial MIP solvers, e.g. GUROBI (Gurobi, 2023), CPLEX (IBM, 2023). In this case, a MIP solver provides the optimality gap (MIP gap) that gauges the quality of the approximate solution, i.e. the last feasible solution (incumbent). This optimality gap is calculated based on the incumbent and the so-called MIP best bound.

In the case of the Multi-Objective MIP (MOMIP) problem, scalarization techniques and MIP solvers can be used to derive Pareto optimal solutions (see, e.g. Miettinen, 1999; Ehrgott, 2005). Examples of applying MIP packages to solve multi-criteria decision problems are shown in, e.g. Ahmadi et al. (2012), Delorme et al. (2014), Eiselt and Marianov (2014), Oke and Siddiqui (2015), Samanlioglu (2013). As a scalarization technique, one can use the Chebyshev scalarization that guarantees the derivation of each (properly) Pareto optimal solution (see, e.g. Kaliszewski, 2006). Other advantages of using this scalarization in the context of decision-making and expressing the decision maker’s preferences are discussed in, e.g. Miroforidis (2021).

In the current work, we say that an instance of the MOMIP problem is large-scale if its Chebyshev scalarization cannot be solved to optimality by a MIP solver within an assumed time limit that is reasonable in the decision-making process. The existence of this limit is justified in solving practical multi-criteria decision-making problems. When there is a time limit on deriving a single Pareto optimal solution, the Chebyshev scalarization of the instance may not be solved to optimality. The decision maker (DM) then obtains the incumbent, i.e. the approximation of the Pareto optimal solution, as well as the MIP gap of the single-objective optimization problem. However, based on this information, the quality of the approximation of a single component (namely its lower and upper bounds) of the Pareto optimal outcome, i.e. the image of the Pareto optimal solution in the objective space cannot be shown to the DM. And it is based on these components that the DM navigates on the Pareto front (set of Pareto optimal outcomes). Fortunately, there is a method to provide the DM with such lower and upper bounds in the literature.

In Kaliszewski and Miroforidis (2019), a general methodology for multi-objective optimization to provide lower and upper bounds on objective function values of a Pareto optimal solution designated by a vector of weights of the Chebyshev scalarization of a multi-objective optimization problem has been proposed. The bounds form the so-called interval representation of the Pareto optimal outcome. The DM can use interval representations instead of (unknown to him/her) Pareto optimal outcomes, to navigate on the Pareto front. To derive them, one needs the so-called lower shells and upper shells whose images in the objective space are finite two-sided approximations of the Pareto front (see, e.g. Kaliszewski and Miroforidis, 2014).

In Kaliszewski and Miroforidis (2022), it has been shown how to provide lower and upper shells to large-scale instances of the MOMIP problem. In that work, lower shells are composed of incumbents to the Chebyshev scalarization of the MOMIP problem derived within the time limit, and upper shells consist of elements that are solutions to the Chebyshev scalarization of a relaxation of the MOMIP problem.

However, there is a lack of an algorithmic method for deriving an upper shell that is necessary to calculate the interval representation of the Pareto optimal outcome designated by a given vector of weights of the Chebyshev scalarizing function. The idea of how to derive such useful upper shells for the MOMIP problem with two objective functions has been shown in our earlier works (Kaliszewski and Miroforidis, 2021) and (Miroforidis, 2021).

In the current work, we combine ideas from works (Kaliszewski and Miroforidis, 2019, 2021, 2022), and (Miroforidis, 2021). For this reason, our work is an incremental one. For the MOMIP problem with up to three objective functions, we propose an algorithmic method of deriving upper shells that can be used to calculate the interval representation of a single Pareto optimal outcome designated by a given vector of weights of the Chebyshev scalarizing function. This opens the way for providing the DM with this representation when there is a time limit for deriving a single Pareto optimal solution. Because of the need to derive the appropriate upper shells, additional time is needed for optimization, but as we show in numerical experiments, this time can be a fraction of the assumed time limit. To illustrate our method, we present results of several numerical experiments with selected large-scale instances of the multi-objective multidimensional 0–1 knapsack problem.

To our best knowledge, the method we propose is the only algorithmic method for determining the interval representation of the Pareto optimal outcome given by weights of the Chebyshev scalarizing function for large-scale instances of the MOMIP problem, assuming the existence of a time budget for optimization.

The main contribution of this article is summarized as follows:

• We propose a generic framework for providing interval representations of Pareto optimal outcomes, designated by weights of the Chebyshev scalarization, of the MOMIP problem when there is a time budget for optimization.
• We propose two algorithms, which are realizations of the framework.
• We demonstrate the operation of these algorithms on computationally demanding large-scale instances of the MOMIP problem up to three objective functions.
• We discuss possible directions for changing the framework to better adapt it to the realities of decision-making and the budgeting of calculations.

The current work is organized as follows. In Section 2, we formulate the MOMIP problem and we recall a method for the derivation of Pareto optimal solutions with the use of the Chebyshev scalarization. In Section 3, we briefly recall the theory of parametric lower and upper bounds. There, we also introduce the concept of the interval representation of the implicit Pareto optimal outcome as well as an indicator measuring its quality. In Section 4, we present two versions of an algorithm for deriving interval representations of implicit Pareto optimal outcomes. In Section 5, we conduct extensive numerical experiments, as well as discuss their results. In Section 6, we show the limitations of the proposed method, as well as discuss how to eliminate them. Section 7 contains some final remarks.

2Background

In this section, we formulate the MOMIP problem, and we recall a method for the derivation of Pareto optimal solutions with the use of the Chebyshev scalarization.

Let X:=Rn1×Zn2, n1+n2=n, n2>0, x∈X denote a solution, X0:={x∈X|gp(x)⩽bp,bp∈R} the set of feasible solutions, where gp:Rn→R, p=1,…,m, m⩾1. The MOMIP problem is defined as follows:

(1)

vmaxf(x)s.t.x∈X0,

where f:Rn→Rk, f=(f1,…,fk), fl:Rn→R, l=1,…,k, k⩾2, are objective functions, and “vmax” is the operator of deriving set N that contains all Pareto optimal solutions in X0. The set Rk is called the objective space. Solution x¯∈X0 is Pareto optimal, if for any x∈X0, fl(x)⩾fl(x¯), l=1,…,k, implies f(x)=f(x¯). If fl(x)⩾fl(x¯), l=1,…,k, and f(x)≠f(x¯), then x dominates x¯ (x¯ is dominated) which is denoted by the relation x≻x¯. We say that element f(x), x∈X0, is the outcome of x. Set f(N) is called the Pareto front.

According to well-established knowledge (Ehrgott, 2005; Kaliszewski, 2006; Kaliszewski et al., 2016; Miettinen, 1999), solution x is Pareto optimal (actually, x is properly Pareto optimal, see, e.g. Ehrgott, 2005; Kaliszewski, 2006; Kaliszewski et al., 2016; Miettinen, 1999) if and only if it solves the Chebyshev scalarization of problem (1), namely

(2)

minmaxlλl(yl∗−fl(x))+ρek(y∗−f(x))s.t.x∈X0,

where weights λl>0, l=1,…,k, ek=(1,1,…,1), yl∗=yˆl+ε, yˆl=maxx∈X0fl(x), yˆl<∞, l=1,…,k, ε>0, and ρ is a positive “sufficiently small” number.

The linearized version of problem (2) is the following.

(3)

minss.t.s⩾λl(yl∗−fl(x))+ρek(y∗−f(x)),l=1,…,k,x∈X0.

In the following, we will assume that Pareto optimal solutions come from solving problem (3) with varying λ=(λ1,…,λk).

Given λ, xPopt(λ) denotes the implicit Pareto optimal solution designated by λ that is a solution which would be derived if problem (2) with λ were solved to optimality. f(xPopt(λ)) denotes the implicit Pareto optimal outcome designated by λ.

3Lower and Upper Bounds on Components of Implicit Pareto Optimal Outcomes

This section contains a brief description of the general theory of lower and upper bounds on components of implicit Pareto optimal outcomes proposed in Kaliszewski and Miroforidis (2019).

To calculate the bounds, one needs two finite sets (that satisfy certain properties) namely a lower shell (SL⊆X0) and upper shell (SU⊆Rn).

Given λ, SL, and SU, the theory provides formulas for calculating lower and upper bounds on fl(xPopt(λ)), l=1,…,k. That is,

(4)

Ll(SL,λ)⩽fl(xPopt(λ))⩽Ul(SU,λ),l=1,…,k.

Formulas for lower bounds Ll(SL,λ) and upper bounds Ul(SU,λ) are shown in Kaliszewski and Miroforidis (2022). In that work, all those elements of the theory of lower and upper bounds that are required to understand the rest of the current work are presented in a synthetic way.

In addition, and of great relevance to the current work, the theory specifies that only elements x∈SU appropriately located with respect to the vector of lower bounds L(SL,λ):=(L1(SL,λ),…,Lk(SL,λ)) can provide upper bounds Ul¯({x},λ)=fl¯(x) for fl¯(xPopt(λ)) for some l¯. This is specified by the following lemma defined in Kaliszewski and Miroforidis (2019).

Lemma 1.

Given lower shell SL and upper shell SU, suppose x∈SU and Ll¯(SL,λ)⩽fl¯(x) for some l¯ and Ll(SL,λ)⩾fl(x) for all l=1,…,k, l≠l¯. Then x provides an upper bound for fl¯(xPopt(λ)), namely fl¯(xPopt(λ))⩽fl¯(x).

Let S¯U⊆SU be a set of elements fulfilling Lemma 1 for some l¯∈{1,…,k}. If S¯U≠∅, then each x∈S¯U can provide an upper bound on fl¯(xPopt(λ)), and Ul¯(SU,λ)=minx∈S¯UUl¯({x},λ). If S¯U=∅, then Ul¯(SU,λ)=yˆl¯. In Section 5, we show how to set Ul¯(SU,λ) when yˆl¯ is not known.

Fig. 1

Components of R(SL,SU,λ): □, ∘ – images of lower shell SL and upper shell SU elements, respectively, in the objective space, ▲ – vector of lower bounds, ■ – vector of upper bounds.

Further on, U(SU,λ):=(U1(SU,λ),…,Uk(SU,λ)) denotes the vector of upper bounds.

3.1The Interval Representation of the Implicit Pareto Optimal Outcome

Given λ, SL, and SU, the interval representation of f(xPopt(λ)) is R(SL,SU,λ)=([L1(SL,λ),U1(SU,λ)],…,[Lk(SL,λ),Uk(SU,λ)]).

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

To gauge the quality of R(SL,SU,λ), we calculate

(5)

GPsub,l(λ):=100×Ul(SU,λ)−Ll(SL,λ)Ul(SU,λ),l=1,…,k.

GPsub(R(SL,SU,λ)):=(GPsub,1(λ),…,GPsub,k(λ)) forms the Pareto suboptimality gap of interval representation R(SL,SU,λ).

4Providing Interval Representations of Implicit Pareto Optimal Outcomes

In this section, we develop a generic framework for providing interval representations of Pareto optimal outcomes, designated by weights of the Chebyshev scalarization, of the MOMIP problem when there is a time limit for optimization.

Given λ, we assume that there is a time limit TL on solving problem (2) by a MIP solver. We also assume that if the MIP solver can not derive the solution to (2), i.e. xPopt(λ), within TL, then it provides incumbent INCλ that is the approximation of xPopt(λ). In this case, our goal is to provide R(SL,SU,λ) calculated on some lower shell SL and on some upper shell SU.

4.1The Derivation of Lower and Upper Shells

As in Kaliszewski and Miroforidis (2022), we will use SL:={INCλ} as a valid lower shell one can use to calculate Ll(SL,λ), l=1,…,k.

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

Lemma 2.

Given λ′, solution x′ to the relaxation of problem (2) with X0′, X0′⊃X0, is not dominated by solution x to problem (2) for any λ.

Lemma 3.

If x is a Pareto optimal solution to the relaxation of problem (1) with X0′, then set {x} is an upper shell to problem (1).

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

The surrogate relaxation of problem (1) with X0′(μ):={x∈X|∑p=1mμpgp(x) ⩽∑p=1mμpbp} instead of X0, 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 SL, based on Lemma 1, element x of upper shell SU is a source of an upper bound on fl¯(xPopt(λ)), l¯∈{1,…,k}, when f(x) is appropriately located with respect to the vector of lower bounds L(SL,λ). In Miroforidis (2021), an idea of how to derive an upper shell that consists of an element useful to calculate upper bounds on fl¯(xPopt(λ)) 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.

Fig. 2

Deriving upper shell SU1, whose element xλ‴ is a source of an upper bound, U1, for f1(xPopt(λ)) with some λ : ∘ – image of upper shell SU1 in the objective space, ▲ – vector of lower bounds.

For k=2, the idea of an algorithm for deriving upper shell SU1 whose some element is a source of an upper bound on f1(xPopt(λ)) is shown in Fig. 2. Let us assume that vector μ is given. At the beginning, SU1:=∅. 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(X0′(μ),λ′). SU1:=SU1∪{xλ′}. Based on Lemma 1, xλ′ is not a source of an upper bound on f1(xPopt(λ)) because f(xλ′) is not appropriately located with respect to the vector of lower bounds L=L(SL,λ). Hence, we CONTINUE.

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

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

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

To obtain an upper bound on f2(xPopt(λ)), we need to derive upper shell SU2. 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 SUl¯ whose element x is a source of an upper bound on fl¯(xPopt(λ)), i.e. Ul¯(SUl¯,λ).

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 SUl¯ checking if its new element fulfills conditions of Lemma 1 to be a valid source for the upper bound on fl¯(xPopt(λ)). 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=1kλl′=1. The exit condition of the “while” loop ensures that λl′>0, l=1,…,k.

If no element of SUl¯ satisfies conditions of Lemma 1, then the algorithm simply returns the upper shell as well as the only available upper bound on fl¯(xPopt(λ)), namely yl¯∗. Otherwise, an upper bound better than yl¯∗ is returned, as well as the upper shell.

4.2Calculating Interval Representations

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)∈Rm. However, we can use the theory of duality for this purpose.

Given μ, λ, let x be the solution to ChebRLX(X0′(μ),λ), 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 X0. 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{minx∈X0′(μ)maxlλl(yl∗−fl(x))+ρek(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 TS 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)||∈Rm, 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 TS, i.e. in pseudocode, it can be used as a function Suboptimal(λ,N,TS), 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 TS.

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(SL,SU,λ), so tighter upper bounds on components of f(xPopt(λ)), than the Chute1 algorithm. However, it is at the cost of increasing the computation time relative to Chute1 by at most TS. We will check it experimentally in the next section.

Fig. 3

The idea of deriving upper shell {x(μ˜)} whose element is a source of a better upper bound on f2(xPopt(λ)) than the element of upper shell {x(μI)}.

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)∈Rm, and μ˜ be the output of the Suboptimal algorithm (with some N and TS) for some λ. x(μI) is the solution to optimization problem (2) with X0 replaced with X0′(μI), and s(μI) is the objective function value of x(μI). x(μ˜) is the solution to optimization problem (2) with X0 replaced with X0′(μ˜), 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 f2(xPopt(λ)). 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 f1(x(μ˜))<f1(x(μI)) AND f2(x(μ˜))<f2(x(μI)), element x(μ˜) is a source of a better upper bound on f2(xPopt(λ)) 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 f1(x(μ˜))<f1(x(μI)) AND f2(x(μ˜))<f2(x(μI)) does not hold. In this case, if f2(x(μ˜))⩾f2(x(μI)), we obtain no better upper bound on f2(xPopt(λ)). On the other hand, if f1(x(μ˜))⩾f1(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 f2(xPopt(λ)), we get a better upper bound on f2(xPopt(λ)).

5Computational Experiments

In this section, we present the results of two experiments where we apply algorithms Chute1 and Chute2 presented in Section 4.2 to selected instances of the Multi-Objective Multidimensional 0–1 Knapsack Problem (MOMKP) with two and three objective functions. The instances are demanding for modern MIP solvers.

5.1Multi-Objective Multidimensional 0–1 Knapsack Problem

For k>1, the MOMKP is formulated in the following way.

(7)

vmaxf1(x):=∑j=1nc1,jxj…fk(x):=∑j=1nck,jxjs.t.x∈X0:={x|∑j=1nap,jxj⩽bp,p=1,…,m,xj∈{0,1},j=1,…,n},

where all ap,j, cp,j are non-negative. In Kaliszewski and Miroforidis (2022), it has been explained why the MOMKP is NP-hard.

5.2Test Instances of the MOMIP Problem

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 1

5.3Experimental Setting

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 yl∗, 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<yl∗. 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∗.

Table 1

Vectors λ and lower bounds for test problems Bi6.1 and Bi9.1.

No.	λ		L(SL,λ) for Bi6.1		L(SL,λ) for Bi9.1
1	0.055	0.945	114253.29	130251.56	104466.45	118482.17
2	0.116	0.884	116707.61	129508.69	107215.43	117756.82
3	0.733	0.267	125690.15	122399.79	116288.83	110899.20
4	0.397	0.603	122075.81	126638.06	112806.06	115033.24
5	0.439	0.561	122514.80	126139.05	113385.01	114671.51

We set TL:=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 Ll(SL,λ),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(SL,λ) (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(SL,λ) (Tables 2 and 3, respectively).

Table 2

Vectors λ and lower bounds L(SL,λ) for test problem Three6.1.

No.	λ			L(SL,λ)
1	0.187	0.770	0.043	118622.54	128616.78	87944.84
2	0.521	0.324	0.155	124156.03	123541.43	115957.65
3	0.067	0.680	0.253	90480.66	127282.50	121460.00
4	0.359	0.295	0.346	120988.48	121527.94	123559.14
5	0.136	0.078	0.786	115623.82	108141.30	129431.77

Table 3

Vectors λ and lower bounds L(SL,λ) for test problem Three9.1.

No.	λ			L(SL,λ)
1	0.351	0.351	0.298	111876.06	111861.18	109288.07
2	0.243	0.143	0.614	110262.24	103915.65	114504.59
3	0.278	0.494	0.228	110549.31	114387.77	107363.36
4	0.179	0.471	0.350	105139.63	113934.39	110819.31
5	0.407	0.014	0.579	112514.84	0.00	113292.80

For all problem instances, for no vector λ, the selected MIP solver derived the solution to problem (2) in the assumed TL=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 GPsub(R(SL,SU,λ)).

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

• U(SU,λ) – components of the vector of upper bounds;
• GAPPsub% – components of GPsub(R(SL,SU,λ));
• |SU| – the number of elements in the upper shell;
• Time SU (s) – time to derive the upper shell (in seconds); for Chute2, the running time of the Suboptimal algorithm is given in parentheses.

In bold, we indicate the improvement of a single component of GPsub(R(SL,SU,λ)) 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 GPsub(R(SL,SU,λ)) when changing the value of parameter γ to a higher value, namely, we consider changes from γ=10 to γ=30 and from γ=30 to γ=50.

5.4Experiment 1 – Deriving Interval Representations with the Chute1 Algorithm

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

Table 4

Chute1, vectors of upper bounds for test problem Bi6.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10		γ=30		γ=50
1	120964.00	131117.00	120466.00	131117.00	120093.00	131117.00
2	121666.00	131117.00	121666.00	131117.00	121441.00	131117.00
3	127790.00	126078.00	127635.00	125842.00	127646.00	125642.00
4	125252.00	128990.00	124924.00	128990.00	124852.00	128990.00
5	125502.00	128779.00	125335.00	128665.00	125307.00	128710.00

Table 5

Chute1, GAPPsub% for test problem Bi6.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10		γ=30		γ=50
1	5.55	0.66	5.16	0.66	4.86	0.66
2	4.08	1.23	4.08	1.23	3.90	1.23
3	1.64	2.92	1.52	2.74	1.53	2.58
4	2.54	1.82	2.28	1.82	2.22	1.82
5	2.38	2.05	2.25	1.96	2.23	2.00

Table 6

Chute1, values of |SU|, and Time SU(s) for test problem Bi6.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	11	32	52	1.90	8.00	11.49
2	11	33	54	2.02	5.97	11.52
3	8	21	34	2.07	4.26	9.12
4	7	19	31	2.38	5.90	9.50
5	7	19	32	2.08	5.63	8.59

Table 7

Chute1, vectors of upper bounds for test problem Bi9.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10		γ=30		γ=50
1	116289.00	119365.00	116289.00	119365.00	116193.00	119365.00
2	117277.00	119365.00	117050.00	119365.00	117057.00	119365.00
3	119379.88	118456.00	119379.88	118456.00	119379.88	118404.00
4	119379.88	119365.00	119295.00	119365.00	119329.00	119365.00
5	119379.88	119365.00	119379.88	119365.00	119379.88	119365.00

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.

Table 8

Chute1, GAPPsub% for test problem Bi9.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10		γ=30		γ=50
1	10.17	0.74	10.17	0.74	10.09	0.74
2	8.58	1.35	8.40	1.35	8.41	1.35
3	2.59	6.38	2.59	6.38	2.59	6.34
4	5.51	3.63	5.44	3.63	5.47	3.63
5	5.02	3.93	5.02	3.93	5.02	3.93

Table 9

Chute1, values of |SU|, and TimeSU(s) for test problem Bi9.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	12	36	59	3.06	11.22	20.44
2	14	40	67	3.91	10.70	17.21
3	17	51	84	3.82	13.05	21.78
4	20	59	99	59	13.79	24.32
5	20	60	100	60	14.84	21.51

Table 10

Chute1, vectors of upper bounds for test problem Three6.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10			γ=30			γ=50
1	128872.00	131116.00	122914.00	128872.00	131116.00	122471.00	128872.00	131116.00	122387.00
2	127143.00	127325.00	124450.00	127069.00	127325.00	123193.00	127105.00	127325.00	122899.00
3	124359.00	131116.00	131738.00	124359.00	131116.00	131738.00	124217.00	131116.00	131738.00
4	125062.00	125714.00	126639.00	124809.00	125327.00	126342.00	124747.00	125273.00	126099.00
5	128872.00	120016.00	131738.00	128872.00	120016.00	131738.00	128872.00	119610.00	131738.00

Table 11

Chute1, GAPPsub% for test problem Three6.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10			γ=30			γ=50
1	7.95	1.91	28.45	7.95	1.91	28.19	7.95	1.91	28.14
2	2.35	2.97	6.82	2.29	2.97	5.87	2.32	2.97	5.65
3	27.24	2.92	7.80	27.24	2.92	7.80	27.16	2.92	7.80
4	3.26	3.33	2.43	3.06	3.03	2.20	3.01	2.99	2.01
5	10.28	9.89	1.75	10.28	9.89	1.75	10.28	9.59	1.75

Table 12

Chute1, values |SU|, and TimeSU(s) for test problem Three6.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	8	14	33	2.69	4.46	10.48
2	11	21	50	3.39	7.80	23.98
3	11	21	49	5.37	11.84	22.38
4	7	12	28	5.46	9.28	24.29
5	11	21	51	3.27	6.54	16.37

Table 13

Chute1, vectors of upper bounds for test problem Three9.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10			γ=30			γ=50
1	119379.88	119365.00	117516.00	119379.8835	119365.00	117398.00	119379.8835	119365.00	117362.00
2	119379.88	119365.00	118122.00	119379.8835	119365.00	118122.00	119379.8835	119365.00	118122.00
3	119379.88	119365.00	118122.00	119379.8835	119365.00	118122.00	119379.8835	119365.00	118122.00
4	119379.88	119365.00	118122.00	119379.8835	119365.00	118122.00	119379.8835	119365.00	118122.00
5	119379.88	119365.00	118122.00	119379.8835	119365.00	118122.00	119379.8835	119365.00	118122.00

Table 14

Chute1, GAPPsub% for test problem Three9.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10			γ=30			γ=50
1	6.29	6.29	7.00	6.29	6.29	6.91	6.29	6.29	6.88
2	7.64	12.94	3.06	7.64	12.94	3.06	7.64	12.94	3.06
3	7.40	4.17	9.11	7.40	4.17	9.11	7.40	4.17	9.11
4	11.93	4.55	6.18	11.93	4.55	6.18	11.93	4.55	6.18
5	5.75	100.00	4.09	5.75	100.00	4.09	5.75	100.00	4.09

Table 15

Chute1, values of |SU|, and TimeSU(s) for test problem Three9.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	28	79	130	50.09	163.04	237.71
2	18	53	86	46.24	324.61	413.15
3	25	69	115	79.80	222.40	387.76
4	22	64	105	36.09	104.40	175.58
5	11	29	49	11.38	50.33	101.04

For instance Bi6.1, when changing γ=10 to γ=30, we observe an improvement in at least one component of GPsub(R(SL,SU,λ)) 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 GPsub(R(SL,SU,λ)) in at least one of its components for three vectors λ. For λ No. 3, we observe a deterioration of the first component of GPsub(R(SL,SU,λ)), and at the same time, an improvement of the second. For λ No. 5, we observe an improvement of the first component of GPsub(R(SL,SU,λ)), 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 GPsub(R(SL,SU,λ)) 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 GPsub(R(SL,SU,λ)) for just two vectors λ. When changing γ=10 to γ=50, we observe an improvement in GPsub(R(SL,SU,λ)) for vectors λ Nos. 1–4 (at least one component improves). For λ No. 5, GPsub(R(SL,SU,λ)) 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 TL=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 GPsub(R(SL,SU,λ)) 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 GPsub(R(SL,SU,λ)) in at least one of its components for three vectors λ. We also observe a deterioration of the first component of GPsub(R(SL,SU,λ)) for λ No. 2, and, at the same time, an improvement on its third component. When changing γ=10 to γ=50, we observe an improvement in GPsub(R(SL,SU,λ)) for all vectors λ (at least one component improves). For instance Three9.1 and vectors λ Nos. 2–5, upper bounds on all components of f(xPopt(λ)) are equal to the corresponding component of y∗. Only for λ No. 1, f3(xPopt(λ))<y3∗, and when changing γ=10 to γ=30, as well as γ=30 to γ=50, there is an improvement only for the third component of GPsub(R(SL,SU,λ)). 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 TL=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 λ.

We checked that for all instances, for no vector λ, the MIP solver derived the optimal solution to problem (2) within time limit TL+TimeSU.

5.5Experiment 2 – Deriving Interval Representations with the Chute2 Algorithm

In this experiment, we check the behaviour of the Chute2 algorithm with N=20 and TS=400. Recall that these parameters are related to the Suboptimal algorithm.

The results for instance Bi6.1 are shown in Tables 16–18, and the results for instance Bi9.1 are shown in Tables 19–21. The results for instance Three6.1 are shown in Tables 22–24, and the results for instance Three9.1 are shown in Tables 25–27.

Table 16

Chute2, upper bounds for test problem Bi6.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10		γ=30		γ=50
1	118423.00	130703.00	116671.00	130703.00	116230.00	130703.00
2	119507.00	130021.00	118226.00	129980.00	117969.00	129946.00
3	126391.00	123927.00	126213.00	123235.00	126179.00	123296.00
4	123079.00	127228.00	122598.00	127102.00	122657.00	127146.00
5	123474.00	126864.00	123273.00	126864.00	123233.00	126766.00

Table 17

Chute2, GAPPsub% for test problem Bi6.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10		γ=30		γ=50
1	3.52	0.35	2.07	0.35	1.70	0.35
2	2.34	0.39	1.28	0.36	1.07	0.34
3	0.55	1.23	0.41	0.68	0.39	0.73
4	0.82	0.46	0.43	0.37	0.47	0.40
5	0.78	0.57	0.62	0.57	0.58	0.49

Table 18

Chute2, |SU|, and TimeSU(s) for test problem Bi6.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	6	16	26	81.16 (78.64)	84.79 (77.97)	91.44 (79.81)
2	4	9	13	182.12 (180.58)	215.77 (211.45)	197.74 (192.33)
3	3	5	8	37.59 (36.74)	41.65 (38.50)	42.82 (37.99)
4	2	3	6	42.39 (41.04)	41.65 (40.42)	44.57 (41.78)
5	2	5	7	44.07 (43.33)	41.28 (39.66)	39.82 (37.83)

Table 19

Chute2, upper bounds for test problem Bi9.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10		γ=30		γ=50
1	110646.00	119365.00	109109.00	119365.00	109313.00	119365.00
2	111592.00	119218.00	111080.00	119124.00	110985.00	119105.00
3	118331.00	114263.00	118172.00	113760.00	118138.00	113655.00
4	115241.00	116795.00	115230.00	116781.00	115121.00	116732.00
5	115443.00	116518.00	115426.00	116271.00	115321.00	116233.00

Table 20

Chute2, GAPPsub% for test problem Bi9.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10		γ=30		γ=50
1	5.58	0.74	4.25	0.74	4.43	0.74
2	3.92	1.23	3.48	1.15	3.40	1.13
3	1.73	2.94	1.59	2.51	1.57	2.42
4	2.11	1.51	2.10	1.50	2.01	1.46
5	1.78	1.58	1.77	1.38	1.68	1.34

Table 21

Chute2, |SU|, and TimeSU(s) for test problem Bi9.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	11	31	52	439.96 (411.14)	515.48 (402.94)	665.55 (404.57)
2	10	27	44	482.23 (407.04)	678.38 (402.56)	765.19 (409.44)
3	8	20	32	443.50 (406.40)	534.97 (410.38)	667.23 (404.69)
4	5	15	23	424.95 (403.06)	465.13 (400.30)	484.92 (402.56)
5	5	13	20	432.90 (400.51)	485.67 (404.07)	512.03 (403.40)

For instance Bi6.1, when changing γ=10 to γ=30, we observe an improvement in at least one component of GPsub(R(SL,SU,λ)) for all vectors λ, although only in three cases (λ Nos. 3–5) two components improve. Yet, when changing γ=30 to γ=50, we observe an improvement of GPsub(R(SL,SU,λ)) in at least one of its components for three vectors λ (Nos. 1, 2, and 5). For λ No. 3, we observe a deterioration of the second component of GPsub(R(SL,SU,λ)), and, at the same time, an improvement on the first one. All components of GPsub(R(SL,SU,λ)) deteriorate for λ No. 4. When changing γ=10 to γ=50, we observe an improvement in GPsub(R(SL,SU,λ)) for all vectors λ (at least one component improves).

Table 22

Chute2, upper bounds for test problem Three6.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10			γ=30			γ=50
1	128872.00	131116.00	122259.00	128872.00	131116.00	120988.00	128872.00	131116.00	120434.00
2	125477.00	125196.00	121464.00	125483.00	125196.00	120230.00	125477.00	125196.00	119034.00
3	122652.00	131116.00	131738.00	122094.00	131116.00	131738.00	122334.00	131116.00	131738.00
4	122671.00	123772.00	125154.00	122242.00	123201.00	124765.00	122404.00	123097.00	124677.00
5	128872.00	119155.00	131738.00	128872.00	118229.00	131738.00	128872.00	117837.00	131738.00

Table 23

Chute2, GAPPsub% for test problem Three6.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10			γ=30			γ=50
1	7.95	1.91	28.07	7.95	1.91	27.31	7.95	1.91	26.98
2	1.05	1.32	4.53	1.06	1.32	3.55	1.05	1.32	2.58
3	26.23	2.92	7.80	25.89	2.92	7.80	26.04	2.92	7.80
4	1.37	1.81	1.27	1.03	1.36	0.97	1.16	1.27	0.90
5	10.28	9.24	1.75	10.28	8.53	1.75	10.28	8.23	1.75

Table 24

Chute2, |SU|, and TimeSU(s) for test problem Three6.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	8	20	31	106.30 (104.24)	133.88 (121.32)	118.75 (102.00)
2	4	11	17	413.24 (408.25)	436.70 (409.23)	429.80 (401.94)
3	10	27	44	50.87 (47.43)	76.35 (61.08)	60.59 (41.18)
4	3	6	10	293.45 (289.56)	380.25 (365.56)	353.60 (319.88)
5	11	30	50	56.90 (51.22)	66.93 (51.58)	79.55 (50.46)

Table 25

Chute2, upper bounds for test problem Three9.1 and γ∈{10,30,50}.

No.	U(SU,λ)
No.	γ=10			γ=30			γ=50
1	115765.00	115804.00	113089.00	115250.00	115546.00	112742.00	115126.00	115335.00	112656.00
2	114092.00	113121.00	118122.00	113842.00	112477.00	117089.00	113861.00	112778.00	118122.00
3	116240.00	117966.00	112426.00	116229.00	117577.00	111977.00	116160.00	117657.00	111578.00
4	112893.00	116633.00	114061.00	112291.00	116382.00	113713.00	112177.00	116438.00	113665.00
5	119379.88	112286.00	118122.00	119379.88	111440.00	118122.00	119379.88	111242.00	118122.00

Table 26

Chute2, GAPPsub% for test problem Three9.1 and γ∈{10,30,50}.

No.	GAPPsub%
No.	γ=10			γ=30			γ=50
1	3.36	3.40	3.36	2.93	3.19	3.06	2.82	3.01	2.99
2	3.36	8.14	3.06	3.14	7.61	2.21	3.16	7.86	3.06
3	4.90	3.03	4.50	4.89	2.71	4.12	4.83	2.78	3.78
4	6.87	2.31	2.84	6.37	2.10	2.54	6.27	2.15	2.50
5	5.75	100.00	4.09	5.75	100.00	4.09	5.75	100.00	4.09

Table 27

Chute2, |SU|, and TimeSU(s) for test problem Three9.1 and γ∈{10,30,50}.

No.	\|SU\|			TimeSU(s)
No.	γ=10	γ=30	γ=50	γ=10	γ=30	γ=50
1	8	20	31	519.63 (420.77)	548.70 (400.11)	654.78 (411.71)
2	12	32	55	412.71 (403.34)	478.80 (409.01)	521.99 (401.65)
3	12	32	52	444.12 (400.92)	527.00 (404.25)	605.74 (405.75)
4	9	22	36	425.02 (400.11)	462.81 (405.86)	515.93 (405.72)
5	5	12	20	325.93 (323.31)	335.62 (324.23)	360.17 (330.95)

For instance Bi9.1, when changing γ=10 to γ=30, we observe an improvement in at least one component of GPsub(R(SL,SU,λ)) for all vectors λ, and for λ Nos. 2–5, both components of GPsub(R(SL,SU,λ)) improve. When changing γ=30 to γ=50, we observe an improvement in both components of GPsub(R(SL,SU,λ)) for all vectors λ but the first one, where the first component of GPsub(R(SL,SU,λ)) deteriorates. When changing γ=10 to γ=50, we observe an improvement in GPsub(R(SL,SU,λ)) for all vectors λ (at least one component improves).

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 instance Bi6.1, average times over all vectors λ to derive the upper shell for γ=10, γ=30, and γ=50 are, respectively, 77.74, 85.03, and 83.28 seconds. These times are small fractions of the assumed time limit TL=1200 seconds. For instance Bi9.1, average times over all vectors λ to derive the upper shell for γ=10, γ=30, and γ=50 are, respectively, 444.71, 535.92, and 618.98 seconds, and they are not small fractions of TL=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 GPsub(R(SL,SU,λ)) in at least one of its components for vectors λ Nos. 1, and 3–5. For λ No. 2, we observe an improvement in the third component of GPsub(R(SL,SU,λ)), as well as a deterioration of the first one. When changing γ=30 to γ=50, we observe an improvement of GPsub(R(SL,SU,λ)) in at least one of its components for vectors λ Nos. 1–3, and 5. For λ No. 4, we observe an improvement in the second and third components of GPsub(R(SL,SU,λ)), as well as a deterioration of the first one. When changing γ=10 to γ=50, we observe an improvement in the GPsub(R(SL,SU,λ)) for all vectors λ (at least one component improves).

For instance Three9.1, when changing γ=10 to γ=30, we observe an improvement in all components of GPsub(R(SL,SU,λ)) for vectors λ Nos. 1–4, and for λ No. 5, GPsub(R(SL,SU,λ)) remains unchanged. When changing γ=30 to γ=50, we observe an improvement in all components of GPsub(R(SL,SU,λ)) for vector λ No. 1. For λ No. 2, we observe a deterioration of all components of GPsub(R(SL,SU,λ)), and for λ Nos. 3–4, we observe an improvement in two components of GPsub(R(SL,SU,λ)), and a deterioration of one of its components. When changing γ=10 to γ=50, we observe an improvement in the GPsub(R(SL,SU,λ)) for all vectors λ (at least one component improves) but the last one.

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 Three9.1, average times over all vectors λ to derive the upper shell for γ=10, γ=30, and γ=50 are, respectively, 425.48, 470.58, and 531.72 seconds, and they are not small fractions of TL=1200 seconds.

We checked that for all instances, for no vector λ, the MIP solver derived the optimal solution to problem (2) within time limit TL+TimeSU.

5.6Comparing Chute2 with Chute1

When comparing Chute2 to Chute1, for all tested instances and all values of the γ parameter, we observe no deterioration of any component of GPsub(R(SL,SU,λ)). We observe the following.

For instance Bi6.1, for all values of γ, all components of GPsub(R(SL,SU,λ)) improve for all vectors λ.

For instance Bi9.1, for γ=10, all components of GPsub(R(SL,SU,λ)) improve for four vectors λ, and one component improves for one vector λ. The same situation occurs for γ=30 and γ=50.

For instance Three6.1, for γ=10, at least one component of GPsub(R(SL,SU,λ)) improves for all vectors λ, and for two ones all components improve. The same situation occurs for γ=30 and γ=50.

For instance Three9.1, for all γ, at least one component of GPsub(R(SL,SU,λ)) improves for four vectors λ. All components of GPsub(R(SL,SU,λ)) improve for γ=10, γ=30, and γ=50, respectively, for three, four, and three vectors λ. For all values of the γ parameter, for one vector λ there is no improvement of any component of GPsub(R(SL,SU,λ)).

Table 28 shows times of deriving upper shells averaged over all vectors λ for both tested algorithms. We observe a significant increase in these times for Chute2 compared to Chute1. It should be recalled here that Chute2 uses the Suboptimal algorithm, for which the stopping condition depends on the assumed for this algorithm time limit TS=400 seconds. For Chute2, the average running time of the Suboptimal algorithm is given in parentheses. It can be seen that in this case a significant fraction of its running time is that of the Subotimal algorithm. In addition, for all γ values, the average running times of Chute2 are larger for Bi9.1 than for Three9.1, which is theoretically a harder problem to solve because it has one more objective function. The implication is that for Bi9.1, for all five lambda vectors, the Subotimal algorithm terminated due to the TS limit, while for Three9.1 – for four lambda vectors. This affected the average times. For details, see Tables 21 and 27.

Table 28

Average times of deriving upper shells for Chute1 and Chute2.

γ	AVG TimeSU(s)	AVG TimeSU(s)
γ	Chute1	Chute2
Bi6.1
10	2.09	77.47 (76.07)
30	5.95	85.03 (81.60)
50	10.04	83.28 (77.95)
Bi9.1
10	4.02	444.71 (405.63)
30	12.72	535.92 (404.05)
50	21.05	618.98 (404.93)
Three6.1
10	4.04	184.15 (180.14)
30	7.98	218.82 (201.75)
50	19.50	208.46 (183.09)
Three9.1
10	44.72	425.48 (389.69)
30	172.96	470.58 (388.69)
50	263.05	531.72 (391.16)

5.7Discussion

For all test instances of the MOMIP problem, with time limits set, algorithm Chute2 determines tighter upper bounds measured with the help of GPsub(R(SL,SU,λ)) than algorithm Chute1 in most cases. Yet, this comes at the expense of a significant increase in the computation time for deriving upper shells. So, we can observe a trade-off between the quality of the interval representation of the implicit Pareto optimal outcome for a given λ and computation time. In both the algorithms, for a given λ, tightness of upper bounds can be controlled by changing values of parameter γ. However, changing the γ value from lower to higher does not always guarantee an improvement in at least one component of GPsub(R(SL,SU,λ)). It may even happen that some of its components deteriorate. However, in all tested instances, when changing from the lowest to the highest value of parameter γ, no deterioration of any component of GPsub(R(SL,SU,λ)) has been recorded for all vectors λ.

The deterioration of some of the components of GPsub(R(SL,SU,λ)) after increasing γ may be due to the fact that increasing the value of γ does not preserve the elements of the set SU obtained for smaller γ, but generates a new, denser set SU, yet different in general. These new SU elements may not be able to generate always better, but in some cases generate even slightly worse vectors of upper bounds than those obtained for smaller γ. During decision-making, one can store all the derived upper shells and use their elements in the Chute algorithm as helpers to determine tighter upper bounds for a given λ when the DM asks for them.

Parameters affecting the operation of algorithms Chute1 and Chute2 (in particular, time limits for optimization, as well as parameter γ) were arbitrarily set for the numerical experiments conducted on the selected test instances. We can not recommend the adopted parameter values (e.g. TL=1200 seconds, TS=400 seconds) for other instances of the MOMIP problem. The values of these parameters might depend on the problem to be solved, the available computational resources and the conditions of the decision-making process itself.

As Chute1 and Chute2 use a MIP solver as a black box, it is difficult to provide their theoretical performance, especially since they can work with any instance of the MOMIP problem that meets the very generic assumptions made in this work. During their operation, multiple instances of the single-objective MIP problem are solved, which are parameterized by λ in the case of Chute1 and (λ,μ) in the case of Chute2. Moreover, Chute2 uses the Suboptimal algorithm as a black box, and it is difficult to predict which termination condition of Suboptimal will occur as it runs for different instances of the single-objective MIP problem parametrized by λ.

The Chute algorithm returns not only the interval representation but also lower and upper shells. Let us assume that for a given set Λ:={λ1,λ2,λ3} the algorithm derives upper shells SU(λ1), SU(λ2), and SU(λ3). SU:=⊕ˆi=13SU(λi) (where “⊕ˆ” is an operator of adding two sets and removing dominating elements) is an upper shell, and SL:=⊕i=13{INCλi} (where “⊕” is an operator of adding two sets and removing dominated elements) is a lower shell. One can use SL and SU to calculate interval representations of implicit Pareto optimal outcomes designated by λ∉Λ. For test problem Bi6.1 and γ:=50, images of the lower and upper shells obtained this way (for five considered vectors λ) are shown in Fig. 4. These images form a finite two-sided approximation of the Pareto front. The approximation does not fully cover the entire Pareto front, as it was derived to determine interval representations of implicit Pareto optimal outcomes designated by just the selected five vectors λ.

Fig. 4

A finite two-sided approximation of the Pareto front: □, ∘ – images of lower shell SL and upper shell SU elements in the objective space, respectively, ∙ – y∗.

We can say that we obtained the two-sided approximation of the Pareto front, shaped by the DM’s preferences expressed with the help of vectors λ.

Although we aim not to derive approximations of the entire Pareto front (as in multi-objective branch and bound, see, e.g. Przybylski and Gandibleux, 2017; Forget et al., 2022) for k=2,3, with a fairly large set of evenly distributed vectors λ, one would imagine the corridor in which the Pareto front is located.

6Limitations of the Chute Algorithm and Its Possible Enhancements

In this section, we discuss the limitations and possible enhancements of the Chute algorithm to better adapt it to the realities of decision-making and the budgeting of calculations.

In the Chute algorithm, we have assumed that for all probing vectors λ′ in the FindUpperShell algorithm, the same vector of multipliers μ is used. The Chute1 version inherently uses a single vector μ. Yet, for the Chute2 version, it is just a heuristic assumption that vector μ, set with the help of the Suboptimal algorithm for a given vector λ in line 4 of the Chute algorithm, provides a tight lower bound on values of the objective function of problem (2) for probing vectors λ′ close to λ. However, this need not be the case, especially for vectors λ′, which are significantly different from λ (i.e. when they indicate a significantly different search direction in the objective space).

However, one can imagine version Chute3 of the Chute algorithm in which the determination of vector μ takes place in the FindUpperShell algorithm for each probing λ′ considered in it (or, e.g. for λ′ not close, in a sense, to a given λ). This, at the same time, would require adopting a reasonable time limit on optimization in the Suboptimal algorithm, as we expect many probing vectors λ′ in the FindUpperShell algorithm. This time limit could be, e.g. a fraction of time TS adopted in Chute2. Since the number of probing vectors λ′ is not a priori known, this time limit would have to be determined by some heuristic rule. It is not desirable that excessive time to determine all vectors μ be a barrier to the applicability of the proposed method.

In a real decision-making process using the Chute algorithm, it is possible to calculate a more adjusted value of parameter γ for a new vector λ based on the properties of the lower and upper shells obtained for previous vectors λ, and with which this algorithm was called. Let us look, for example, at Table 18. For λ1=(0.055,0.945), |SU|=26, TimeSU=91.44 seconds, but for λ5=(0.439,0.561), |SU|=7, TimeSU=39.82 seconds. To have the time of deriving an upper shell for some λ′ close to λ1 comparable to the time for λ5, it could be possible to lower the value of parameter γ from 50 to, e.g. 20. Such an overarching mechanism (with a set of rules based on statistics collected during the decision-making process) for controlling the behaviour of the Chute algorithm could be useful when computation time is an important factor.

In the proposed method of deriving upper shells (algorithm FindUpperShell), there is no single parameter to limit optimization time for getting theirs. Yet, such a time limit can be incorporated relatively easily as an additional stop condition in FindUpperShell. More generally, it could be even desirable to introduce in the Chute algorithm a time limit for determining the interval representation of the implicit Pareto optimal outcome for a given vector λ. The DM would give, for example, in addition to λ and TL, time limit TI (e.g. TL=1200 seconds, TI=500 seconds). Then the Chute algorithm would have time limit TIk to derive an upper shell for calculating a single component of the interval representation of implicit Pareto optimal outcome designated by λ. Determination of suboptimal vectors μ in Chute2 and Chute3 would, of course, have to be within some fraction of TI.

Based on the above alone, one can imagine many schemes for budgeting calculations, leading to providing interval representations in a decision-making system based on the Chute algorithm.

In our approach, to find elements of the upper shell we solve to optimality the Chebyshev scalarization of the surrogate relaxation of the MOMIP problem. For instances of the MOMIP problem with a large number of constraints (e.g. 1000), even with a suboptimal vector of multipliers μ provided by the Suboptimal algorithm (that is, with a single constraint that mimics the original set of constraints of the MOMIP problem), the FindUpperShell procedure may not derive elements x of the upper shell that fl(x)<yl∗, l=1,…,k. In this case, the upper bounds on components of Pareto optimal outcome designated by λ are not better than components of yl∗. That is, images of elements of the upper shell in the objective space are very far from the Pareto front of the MOMIP problem, and do not provide better upper bounds than the components of y∗.

To find (sub)optimal values of multipliers μp, other algorithms can be used (see, e.g. Sikorski, 1986). To find elements of upper shells, sophisticated combined relaxation techniques for MIP problems, e.g. Lagrangean/surrogate heuristics (see Narciso and Lorena, 1999) can also be applied. In the current work, we consider the most general formulation of the MOMIP problem, but to find those elements, problem-specific techniques may help. The disadvantage of the proposed generic scheme Chute is that it does not take into account the specifics of a given instance of the MOMIP problem. However, by showing its Chute2 modification and pointing to the Chute3 option, it has been shown how this scheme can be modified.

Within the generic framework presented, other methods of deriving upper shells in the FindUpperShell procedure can also be applied, e.g. a method shown in Miroforidis (2021).

7Final Remarks

It has been shown how to algorithmically derive lower and upper shells to the MOMIP problem (for any k>1) to get the interval representation of the Pareto optimal outcome designated by vector λ. On selected examples, it has been shown that with the help of the proposed method, one can find such interval representations for randomly selected vectors λ where there is a time limit for a MIP solver on deriving a single Pareto optimal solution.

We conducted some preliminary experiments with the Chute algorithm on instances of the MOMKP with four objective functions. However, due to the mechanism adopted in the FindUpperShell algorithm for changing the probing λ vectors, the results achieved were not satisfactory.

In our future work, we want to improve the method of populating upper shells (in quest of finding their elements that can provide upper bounds) by changing the scheme of probing the objective space. We want it to determine upper shells with the desired properties for four and more objective functions. We also want to apply the presented generic approach to other instances of the MOMIP problem, especially ones connected to real-life problems. This would help verify the practicality of the proposed general method and identify those elements that could be tailored for specific instances of this problem. Possible modifications to the proposed method are indicated in Section 6. These are also worth considering in further work.

Notes

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

Acknowledgements

We thank two anonymous reviewers for their helpful comments.

References

1	Ahmadi, A., Aghaei, J., Shayanfar, H.A., Rabiee, A. ((2012) ). Mixed integer programming of multiobjective hydro-thermal self scheduling. Applied Soft Computing, 12: , 2137–2146. https://doi.org/10.1016/j.asoc.2012.03.020.
2	Delorme, X., Battaïa, O., Dolgui, A. ((2014) ). Multi-objective approaches for design of assembly lines. In: Benyoucef, L., Hennet, C., J., Tiwari, M. (Eds.), Applications of Multi-Criteria and Game Theory Approaches: Manufacturing and Logistics. Springer, London, pp. 31–56. https://doi.org/10.1007/978-1-4471-5295-8_2.
3	Dyer, M.E. ((1980) ). Calculating Surrogate Constraints. Mathematical Programming, 19: , 255–278. https://doi.org/10.1007/BF01581647.
4	Ehrgott, M. ((2005) ). Multicriteria Optimization. Springer, Berlin, Heidelberg. https://doi.org/10.1007/3-540-27659-9.
5	Eiselt, H.A., Marianov, V. ((2014) ). A bi-objective model for the location of landfills for municipal solid waste. European Journal of Operational Research, 235: 1, 187–194. https://doi.org/10.1016/j.ejor.2013.10.005.
6	Forget, N., Gadegaard, S.L., Klamroth, K., Nielsen, L.R., Przybylski, A. ((2022) ). Branch-and-bound and objective branching with three or more objectives. Computers & Operations Research, 148: , 106012. https://doi.org/10.1016/j.cor.2022.106012.
7	Glover, F. ((1965) ). A multiphase-dual algorithm for the zero-one integer programming problem. Operations Research, 13: (6), 879–919. https://doi.org/10.1287/opre.13.6.879.
8	Glover, F. ((1968) ). Surrogate constraints. Operations Research, 16: 4, 741–749. https://doi.org/10.1287/opre.16.4.741.
9	Gurobi (2023). GUROBI. Last accessed March 31, 2023. https://www.gurobi.com/.
10	IBM (2023). CPLEX. Last accessed March 31, 2023. https://www.ibm.com/products/ilog-cplex-optimization-studio.
11	Kaliszewski, I. ((2006) ). Soft Computing for Complex Multiple Criteria Decision Making. Springer, New York. https://doi.org/10.1007/0-387-30177-1.
12	Kaliszewski, I., Miroforidis, J. ((2014) ). Two-sided Pareto front approximations. Journal of Optimization Theory and Its Applications, 162: , 845–855. https://doi.org/10.1007/s10957-013-0498-y.
13	Kaliszewski, I., Miroforidis, J. ((2019) ). Lower and upper bounds for the general multiobjective optimization problem. AIP Conference Proceedings, 2070: , 020038. https://doi.org/10.1063/1.5090005.
14	Kaliszewski, I., Miroforidis, J. ((2021) ). Cooperative multiobjective optimization with bounds on objective functions. Journal of Global Optimization, 79: , 369–385. https://doi.org/10.1007/s10898-020-00946-4.
15	Kaliszewski, I., Miroforidis, J. ((2022) ). Probing the Pareto front of a large-scale multiobjective problem with a MIP solver. Operational Research, 22: , 5617–5673. https://doi.org/10.1007/s12351-022-00708-y.
16	Kaliszewski, I., Miroforidis, J., Podkopaev, D. ((2016) ). Multiple Criteria Decision Making by Multiobjective Optimization – A Toolbox. Springer Cham. https://doi.org/10.1007/978-3-319-32756-3.
17	Miettinen, K.M. ((1999) ). Nonlinear Multiobjective Optimization. Kluwer Academic Publishers. https://doi.org/10.1007/978-1-4615-5563-6.
18	Miroforidis, J. ((2021) ). Bounds on efficient outcomes for large-scale cardinality-constrained Markowitz problems. Journal of Global Optimization, 80: , 617–634. https://doi.org/10.1007/s10898-021-01022-1.
19	Narciso, M.G., Lorena, L.A.N. ((1999) ). Lagrangean/surrogate relaxation for generalized assignment problems. European Journal of Operational Research, 114: (1), 165–177. https://doi.org/10.1016/S0377-2217(98)00038-1.
20	Oke, O., Siddiqui, S. ((2015) ). Efficient automated schematic map drawing using multiobjective mixed integer programming. Computers & Operations Research, 61: , 1–17. https://doi.org/10.1016/j.cor.2015.02.010.
21	Przybylski, A., Gandibleux, X. ((2017) ). Multi-objective branch and bound. European Journal of Operational Research, 260: (3), 856–872. https://doi.org/10.1016/j.ejor.2017.01.032.
22	Samanlioglu, F. ((2013) ). A multi-objective mathematical model for the industrial hazardous waste location-routing problem. European Journal of Operational Research, 226: (2), 332–340. https://doi.org/10.1016/j.ejor.2012.11.019.
23	Sikorski, J. ((1986) ). Quasi–Subgradient algorithms for calculating surrogate constraints. In: Malanowski, K., Mizukami, K. (Eds.), Analysis and Algorithms of Optimization Problems, Lecture Notes in Control and Information Sciences, Vol. 82: . Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 202–236. https://doi.org/10.1007/BFb0007163.
24	Smith, N.A., Tromble, R.W. ((2004) ). Sampling Uniformly from the Unit Simplex. Department of Computer Science, Center for Language and Speech Recognition, Johns Hopkins University.