1 Introduction
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.
2 Background
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:={\mathbb{R}^{{n_{1}}}}\times {\mathbb{Z}^{{n_{2}}}}$, ${n_{1}}+{n_{2}}=n$, ${n_{2}}\gt 0$, $x\in X$ denote a solution, ${X_{0}}:=\{x\in X\hspace{0.1667em}|\hspace{0.1667em}{g_{p}}(x)\leqslant {b_{p}},\hspace{2.5pt}{b_{p}}\in \mathbb{R}\}$ the set of feasible solutions, where ${g_{p}}:{\mathbb{R}^{n}}\to \mathbb{R}$, $p=1,\dots ,m$, $m\geqslant 1$. The MOMIP problem is defined as follows:
where $f:{\mathbb{R}^{n}}\to {\mathbb{R}^{k}}$, $f=({f_{1}},\dots ,{f_{k}})$, ${f_{l}}:{\mathbb{R}^{n}}\to \mathbb{R}$, $l=1,\dots ,k$, $k\geqslant 2$, are objective functions, and “vmax” is the operator of deriving set N that contains all Pareto optimal solutions in ${X_{0}}$. The set ${\mathbb{R}^{k}}$ is called the objective space. Solution $\bar{x}\in {X_{0}}$ is Pareto optimal, if for any $x\in {X_{0}}$, ${f_{l}}(x)\geqslant {f_{l}}(\bar{x})$, $l=1,\dots ,k$, implies $f(x)=f(\bar{x})$. If ${f_{l}}(x)\geqslant {f_{l}}(\bar{x})$, $l=1,\dots ,k$, and $f(x)\ne f(\bar{x})$, then x dominates $\bar{x}$ ($\bar{x}$ is dominated) which is denoted by the relation $x\succ \bar{x}$. We say that element $f(x)$, $x\in {X_{0}}$, is the outcome of x. Set $f(N)$ is called the Pareto front.
(1)
\[ \begin{array}{l@{\hskip4.0pt}l}\text{vmax}& f(x)\\ {} \text{s.t.}& x\in {X_{0}},\end{array}\]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
where weights ${\lambda _{l}}\gt 0$, $l=1,\dots ,k$, ${e^{k}}=(1,1,\dots ,1)$, ${y_{l}^{\ast }}={\hat{y}_{l}}+\varepsilon $, ${\hat{y}_{l}}={\max _{x\in {X_{0}}}}{f_{l}}(x)$, ${\hat{y}_{l}}\lt \infty $, $l=1,\dots ,k$, $\varepsilon \gt 0$, and ρ is a positive “sufficiently small” number.
(2)
\[ \begin{array}{l@{\hskip4.0pt}l}\min & {\max _{l}}{\lambda _{l}}\big({y_{l}^{\ast }}-{f_{l}}(x)\big)+\rho {e^{k}}\big({y^{\ast }}-f(x)\big)\\ {} \text{s.t.}& x\in {X_{0}},\end{array}\]The linearized version of problem (2) is the following.
In the following, we will assume that Pareto optimal solutions come from solving problem (3) with varying $\lambda =({\lambda _{1}},\dots ,{\lambda _{k}})$.
(3)
\[ \begin{array}{l@{\hskip4.0pt}l}\min & s\\ {} \text{s.t.}& s\geqslant {\lambda _{l}}({y_{l}^{\ast }}-{f_{l}}(x))+\rho {e^{k}}({y^{\ast }}-f(x)),\hspace{2.5pt}l=1,\dots ,k,\\ {} & x\in {X_{0}}.\end{array}\]Given λ, ${x^{{P_{opt}}}}(\lambda )$ 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({x^{{P_{opt}}}}(\lambda ))$ denotes the implicit Pareto optimal outcome designated by λ.
3 Lower 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 (${S_{L}}\subseteq {X_{0}}$) and upper shell (${S_{U}}\subseteq {\mathbb{R}^{n}}$).
Given λ, ${S_{L}}$, and ${S_{U}}$, the theory provides formulas for calculating lower and upper bounds on ${f_{l}}({x^{{P_{opt}}}}(\lambda ))$, $l=1,\dots ,k$. That is,
Formulas for lower bounds ${L_{l}}({S_{L}},\lambda )$ and upper bounds ${U_{l}}({S_{U}},\lambda )$ 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.
(4)
\[ {L_{l}}({S_{L}},\lambda )\leqslant {f_{l}}\big({x^{{P_{\textit{opt}}}}}(\lambda )\big)\leqslant {U_{l}}({S_{U}},\lambda ),\hspace{1em}l=1,\dots ,k.\]In addition, and of great relevance to the current work, the theory specifies that only elements $x\in {S_{U}}$ appropriately located with respect to the vector of lower bounds $L({S_{L}},\lambda ):=({L_{1}}({S_{L}},\lambda ),\dots ,{L_{k}}({S_{L}},\lambda ))$ can provide upper bounds ${U_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}(\{x\},\lambda )={f_{\bar{l}}}(x)$ for ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ for some $\bar{l}$. This is specified by the following lemma defined in Kaliszewski and Miroforidis (2019).
Lemma 1.
Given lower shell ${S_{L}}$ and upper shell ${S_{U}}$, suppose $x\in {S_{U}}$ and ${L_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({S_{L}},\lambda )\leqslant {f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}(x)$ for some $\bar{l}$ and ${L_{l}}({S_{L}},\lambda )\geqslant {f_{l}}(x)$ for all $l=1,\dots ,k$, $l\ne \bar{l}$. Then x provides an upper bound for ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, namely ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))\leqslant {f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}(x)$.
Let ${\bar{S}_{U}}\subseteq {S_{U}}$ be a set of elements fulfilling Lemma 1 for some $\bar{l}\in \{1,\dots ,k\}$. If ${\bar{S}_{U}}\ne \varnothing $, then each $x\in {\bar{S}_{U}}$ can provide an upper bound on ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, and ${U_{\bar{l}}}({S_{U}},\lambda )={\min _{x\in {\bar{S}_{U}}}}{U_{\bar{l}}}(\{x\},\lambda )$. If ${\bar{S}_{U}}=\varnothing $, then ${U_{\bar{l}}}({S_{U}},\lambda )={\hat{y}_{\bar{l}}}$. In Section 5, we show how to set ${U_{\bar{l}}}({S_{U}},\lambda )$ when ${\hat{y}_{\bar{l}}}$ is not known.
Fig. 1
Components of $R({S_{L}},{S_{U}},\lambda )$: □, ∘ – images of lower shell ${S_{L}}$ and upper shell ${S_{U}}$ elements, respectively, in the objective space, ▲ – vector of lower bounds, ■ – vector of upper bounds.
Further on, $U({S_{U}},\lambda ):=({U_{1}}({S_{U}},\lambda ),\dots ,{U_{k}}({S_{U}},\lambda ))$ denotes the vector of upper bounds.
3.1 The Interval Representation of the Implicit Pareto Optimal Outcome
Given λ, ${S_{L}}$, and ${S_{U}}$, the interval representation of $f({x^{{P_{\textit{opt}}}}}(\lambda ))$ is $R({S_{L}},{S_{U}},\lambda )=([{L_{1}}({S_{L}},\lambda ),{U_{1}}({S_{U}},\lambda )],\dots ,[{L_{k}}({S_{L}},\lambda ),{U_{k}}({S_{U}},\lambda )])$.
For $k=2$, components of the interval representation of $f({x^{{P_{\textit{opt}}}}}(\lambda ))$, 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}},\lambda )$, we calculate
${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda )):=({G_{{P_{sub}},1}}(\lambda ),\dots ,{G_{{P_{sub}},k}}(\lambda ))$ forms the Pareto suboptimality gap of interval representation $R({S_{L}},{S_{U}},\lambda )$.
4 Providing 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 ${T^{L}}$ 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. ${x^{{P_{\textit{opt}}}}}(\lambda )$, within ${T^{L}}$, then it provides incumbent ${INC^{\lambda }}$ that is the approximation of ${x^{{P_{\textit{opt}}}}}(\lambda )$. In this case, our goal is to provide $R({S_{L}},{S_{U}},\lambda )$ calculated on some lower shell ${S_{L}}$ and on some upper shell ${S_{U}}$.
4.1 The Derivation of Lower and Upper Shells
As in Kaliszewski and Miroforidis (2022), we will use ${S_{L}}:=\{{INC^{\lambda }}\}$ as a valid lower shell one can use to calculate ${L_{l}}({S_{L}},\lambda )$, $l=1,\dots ,k$.
To populate upper shell ${S_{U}}$, the following two lemmas defined in Kaliszewski and Miroforidis (2022) can be used.
Given ${X^{\prime }_{0}}\supset {X_{0}}$ and λ, let $\text{ChebRLX}({X^{\prime }_{0}},\lambda )$ denote the Chebyshev scalarization (problem (2)) of some relaxation of the MOMIP problem (problem (1)) with feasible set ${X^{\prime }_{0}}$ for some λ. Based on Lemmas 2 and 3, one can derive a single-element upper shell by solving $\text{ChebRLX}({X^{\prime }_{0}},\lambda )$. Given ${X^{\prime }_{0}}\supset {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^{\prime }_{0}}(\mu ):=\{x\in X\hspace{0.1667em}|\hspace{0.1667em}{\textstyle\sum _{p=1}^{m}}{\mu _{p}}{g_{p}}(x)$ $\leqslant {\textstyle\sum _{p=1}^{m}}{\mu _{p}}{b_{p}}\}$ instead of ${X_{0}}$, where ${\mu _{p}}\geqslant 0$, $p=1,\dots ,m$, $\mu \ne 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_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, $\bar{l}\in \{1,\dots ,k\}$, when $f(x)$ is appropriately located with respect to the vector of lower bounds $L({S_{L}},\lambda )$. In Miroforidis (2021), an idea of how to derive an upper shell that consists of an element useful to calculate upper bounds on ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ 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 ${S_{U}^{1}}$, whose element ${x^{{\lambda ^{\prime\prime\prime }}}}$ is a source of an upper bound, ${U_{1}}$, for ${f_{1}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ with some λ : ∘ – image of upper shell ${S_{U}^{1}}$ in the objective space, ▲ – vector of lower bounds.
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_{\textit{opt}}}}}(\lambda ))$ is shown in Fig. 2. Let us assume that vector μ is given. At the beginning, ${S_{U}^{1}}:=\varnothing $. Below, we describe three iterations of this example.
Iteration 1. We set the first probing vector ${\lambda ^{\prime }}:=({\lambda _{1}}+\delta ,{\lambda _{2}}-\delta )$, ${\lambda ^{\prime }_{2}}\gt 0$, for some $\delta \gt 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^{{\lambda ^{\prime }}}}$ be the solution to $\text{ChebRLX}({X^{\prime }_{0}}(\mu ),{\lambda ^{\prime }})$. ${S_{U}^{1}}:={S_{U}^{1}}\cup \{{x^{{\lambda ^{\prime }}}}\}$. Based on Lemma 1, ${x^{{\lambda ^{\prime }}}}$ is not a source of an upper bound on ${f_{1}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ because $f({x^{{\lambda ^{\prime }}}})$ is not appropriately located with respect to the vector of lower bounds $L=L({S_{L}},\lambda )$. Hence, we CONTINUE.
Iteration 2. We set ${\lambda ^{\prime\prime }}:=({\lambda ^{\prime }_{1}}+\delta ,{\lambda ^{\prime }_{2}}-\delta )$, ${\lambda ^{\prime\prime }_{2}}\gt 0$. Let ${x^{{\lambda ^{\prime\prime }}}}$ be the solution to $\text{ChebRLX}({X^{\prime }_{0}}(\mu ),{\lambda ^{\prime\prime }})$. ${S_{U}^{1}}:={S_{U}^{1}}\cup \{{x^{{\lambda ^{\prime\prime }}}}\}$. Based on Lemma 1, ${x^{{\lambda ^{\prime\prime }}}}$ is not a source of an upper bound on ${f_{1}}({x^{{P_{\textit{opt}}}}}(\lambda ))$. Hence, we CONTINUE.
Iteration 3. We set ${\lambda ^{\prime\prime\prime }}:=({\lambda ^{\prime\prime }_{1}}+\delta ,{\lambda ^{\prime\prime }_{2}}-\delta )$, ${\lambda ^{\prime\prime\prime }_{2}}\gt 0$. Let ${x^{{\lambda ^{\prime\prime\prime }}}}$ be the solution to $\text{ChebRLX}({X^{\prime }_{0}}(\mu ),{\lambda ^{\prime\prime\prime }})$. ${S_{U}^{1}}:={S_{U}^{1}}\cup \{{x^{{\lambda ^{\prime\prime\prime }}}}\}$. Based on Lemma 1, ${x^{{\lambda ^{\prime\prime\prime }}}}$ is a source of an upper bound on ${f_{1}}({x^{{P_{\textit{opt}}}}}(\lambda ))$. So, ${U_{1}}=U({S_{U}^{1}},\lambda )={f_{1}}({x^{{\lambda ^{\prime\prime\prime }}}})$. Hence, we STOP.
As elements ${x^{{\lambda ^{\prime }}}}$, ${x^{{\lambda ^{\prime\prime }}}}$, and ${x^{{\lambda ^{\prime\prime\prime }}}}$ are Pareto optimal solutions to the relaxation of the MOMIP problem with ${X^{\prime }_{0}}(\mu )$, ${S_{U}^{1}}$ is a valid upper shell.
To obtain an upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, we need to derive upper shell ${S_{U}^{2}}$. To do this, in the first iteration, we set ${\lambda ^{\prime }}:=({\lambda _{1}}-\delta ,{\lambda _{2}}+\delta )$, ${\lambda ^{\prime }_{1}}\gt 0$, and proceed in the same way.
Given $\bar{l}\in \{1,\dots ,k\}$, the FindUpperShell algorithm tries to derive upper shell ${S_{U}^{\bar{l}}}$ whose element x is a source of an upper bound on ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, i.e. ${U_{\bar{l}}}({S_{U}^{\bar{l}}},\lambda )$.
In Line 2, we set step size δ that is used to modify components of consecutive probing vectors ${\lambda ^{\prime }}$. Parameter $\gamma \gt 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}^{\bar{l}}}$ checking if its new element fulfills conditions of Lemma 1 to be a valid source for the upper bound on ${f_{\bar{l}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$. The algorithm stops when ${\lambda ^{\prime }_{\bar{l}}}\geqslant 1$ OR some component of ${\lambda ^{\prime }}$ is negative OR an element that fulfills conditions of Lemma 1 is found. Lines 6 and 9 guarantee that ${\textstyle\sum _{l=1}^{k}}{\lambda ^{\prime }_{l}}=1$. The exit condition of the “while” loop ensures that ${\lambda ^{\prime }_{l}}\gt 0$, $l=1,\dots ,k$.
If no element of ${S_{U}^{\bar{l}}}$ satisfies conditions of Lemma 1, then the algorithm simply returns the upper shell as well as the only available upper bound on ${f_{\hspace{0.1667em}\bar{l}\hspace{0.1667em}}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, namely ${y_{\bar{l}}^{\ast }}$. Otherwise, an upper bound better than ${y_{\bar{l}}^{\ast }}$ is returned, as well as the upper shell.
4.2 Calculating 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 ${\mu _{p}}:=1$, $p=1,\dots ,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,\dots ,1)\in {\mathbb{R}^{m}}$. However, we can use the theory of duality for this purpose.
Given μ, λ, let x be the solution to $\text{ChebRLX}({X^{\prime }_{0}}(\mu ),\lambda )$, and $s(\mu )$ 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\geqslant s(\mu )$. Hence, for a given λ and μ, $s(\mu )$ is a lower bound on values of s.
Given λ, the best (highest) lower bound ${s^{\ast }}$ on values of s is the objective function value of the solution ${\mu ^{\ast }}$ to the following surrogate dual problem
that is connected to the Chebyshev scalarization (problem (2)). Solving (6) to optimality can be time-consuming. Yet, a suboptimal vector of multipliers $\tilde{\mu }$ can be determined instead of ${\mu ^{\ast }}$. 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.
(6)
\[ \underset{\mu \geqslant 0,\mu \ne 0}{\sup }\Big\{\underset{x\in {X^{\prime }_{0}}(\mu )}{\min }\underset{l}{\max }{\lambda _{l}}\big({y_{l}^{\ast }}-{f_{l}}(x)\big)+\rho {e^{k}}\big({y^{\ast }}-f(x)\big)\Big\}\]“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 $\frac{(1,\dots ,1)}{||(1,\dots ,1)||}\in {\mathbb{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 $\text{Suboptimal}(\lambda ,N,{T^{S}})$, which returns a suboptimal vector of multipliers $\tilde{\mu }$.
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 ${\lambda ^{\prime }}$, that are close to λ, the Chute2 algorithm is able to find a better $R({S_{L}},{S_{U}},\lambda )$, so tighter upper bounds on components of $f({x^{{P_{\textit{opt}}}}}(\lambda ))$, 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.
Fig. 3
The idea of deriving upper shell $\{x(\tilde{\mu })\}$ whose element is a source of a better upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ than the element of upper shell $\{x({\mu ^{I}})\}$.
The idea (for $k=2$ and $\rho =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 ${\mu ^{I}}:=(1,\dots ,1)\in {\mathbb{R}^{m}}$, and $\tilde{\mu }$ be the output of the Suboptimal algorithm (with some N and ${T^{S}}$) for some λ. $x({\mu ^{I}})$ is the solution to optimization problem (2) with ${X_{0}}$ replaced with ${X^{\prime }_{0}}({\mu ^{I}})$, and $s({\mu ^{I}})$ is the objective function value of $x({\mu ^{I}})$. $x(\tilde{\mu })$ is the solution to optimization problem (2) with ${X_{0}}$ replaced with ${X^{\prime }_{0}}(\tilde{\mu })$, and $s(\tilde{\mu })$ is the objective function value of $x(\tilde{\mu })$. Vector of lower bounds L is marked with a triangle. Both elements $f(x(\tilde{\mu }))$ and $f(x({\mu ^{I}}))$ are appropriately located with regards to L (see Lemma 1) to be sources for an upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$. As $s(\tilde{\mu })\gt s({\mu ^{I}})$ (contours of the Chebyshev metric for both values $s(\tilde{\mu })$ and $s({\mu ^{I}})$ are shown by solid thin lines) as well as ${f_{1}}(x(\tilde{\mu }))\lt {f_{1}}(x({\mu ^{I}}))$ AND ${f_{2}}(x(\tilde{\mu }))\lt {f_{2}}(x({\mu ^{I}}))$, element $x(\tilde{\mu })$ is a source of a better upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$ than element $x({\mu ^{I}})$, as upper bounds are calculated with the use of components of the upper shell elements. In Fig. 3, $f(x(\tilde{\mu }))$ is closer to the (unknown) Pareto front (represented by the solid curve) than element $f(x({\mu ^{I}}))$. It could happen that condition ${f_{1}}(x(\tilde{\mu }))\lt {f_{1}}(x({\mu ^{I}}))$ AND ${f_{2}}(x(\tilde{\mu }))\lt {f_{2}}(x({\mu ^{I}}))$ does not hold. In this case, if ${f_{2}}(x(\tilde{\mu }))\geqslant {f_{2}}(x({\mu ^{I}}))$, we obtain no better upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$. On the other hand, if ${f_{1}}(x(\tilde{\mu }))\geqslant {f_{1}}(x({\mu ^{I}}))$ and still $f(x(\tilde{\mu }))$ is appropriately located with regards to L (see Lemma 1) to be a source for an upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$, we get a better upper bound on ${f_{2}}({x^{{P_{\textit{opt}}}}}(\lambda ))$.
5 Computational 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.1 Multi-Objective Multidimensional 0–1 Knapsack Problem
For $k\gt 1$, the MOMKP is formulated in the following way.
where all ${a_{p,j}}\hspace{0.1667em}$, ${c_{p,j}}$ are non-negative. In Kaliszewski and Miroforidis (2022), it has been explained why the MOMKP is $\mathcal{N}\mathcal{P}$-hard.
(7)
\[ \begin{array}{l@{\hskip4.0pt}l}\text{vmax}& \left\{\begin{array}{l}{f_{1}}(x):={\textstyle\textstyle\sum _{j=1}^{n}}{c_{1,j}}{x_{j}}\\ {} \dots \\ {} {f_{k}}(x):={\textstyle\textstyle\sum _{j=1}^{n}}{c_{k,j}}{x_{j}}\end{array}\right.\\ {} \text{s.t.}& x\in {X_{0}}:=\Bigg\{x\hspace{0.1667em}\Big|\hspace{0.1667em}{\textstyle\textstyle\sum _{j=1}^{n}}{a_{p,j}}{x_{j}}\leqslant {b_{p}},\hspace{2.5pt}p=1,\dots ,m,\\ {} & {x_{j}}\in \{0,1\},\hspace{2.5pt}j=1,\dots ,n\Bigg\},\end{array}\]5.2 Test 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
5.3 Experimental 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 $\hat{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}^{\ast }}$, $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^{\ast }}:=(128872,131116,131738)$, and for Three9.1, ${y^{\ast }}:=(119379.88,119365,118122)$. Obviously, ${\hat{y}_{l}}\lt {y_{l}^{\ast }}$. Hence, such ${y^{\ast }}$ approximating $\hat{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^{\ast }}$.
Table 1
Vectors λ and lower bounds for test problems Bi6.1 and Bi9.1.
| No. | λ | $L({S_{L}},\lambda )$ for Bi6.1 | $L({S_{L}},\lambda )$ 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 ${T^{L}}:=1200$ seconds, $\rho :=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}},\lambda ),l=1,\dots ,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}},\lambda )$ (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}},\lambda )$ (Tables 2 and 3, respectively).
Table 2
Vectors λ and lower bounds $L({S_{L}},\lambda )$ for test problem Three6.1.
| No. | λ | $L({S_{L}},\lambda )$ | ||||
| 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({S_{L}},\lambda )$ for test problem Three9.1.
| No. | λ | $L({S_{L}},\lambda )$ | ||||
| 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 ${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 $\gamma :=10,30,50$. Given λ, we check the impact of parameter γ on components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$.
In the tables with results of the experiments, the meaning of the columns is as follows.
-
• $U({S_{U}},\lambda )$ – components of the vector of upper bounds;
-
• ${GAP_{{P_{sub}}}}\% $ – components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$;
-
• $|{S_{U}}|$ – the number of elements in the upper shell;
-
• Time ${S_{U}}$ (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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ when changing the value of parameter γ to a higher value, namely, we consider changes from $\gamma =10$ to $\gamma =30$ and from $\gamma =30$ to $\gamma =50$. By underscore, we indicate the deterioration of a single component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ when changing the value of parameter γ to a higher value, namely, we consider changes from $\gamma =10$ to $\gamma =30$ and from $\gamma =30$ to $\gamma =50$.
5.4 Experiment 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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Bi6.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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 $|{S_{U}}|$, and Time ${S_{U}}(s)$ for test problem Bi6.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Bi9.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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 $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Bi9.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Three6.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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 $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Three6.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Three9.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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 $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Three9.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma =10$ to $\gamma =30$, we observe an improvement in at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for four vectors λ, although only in two cases (λ Nos. 3 and 5) two components improve. Yet, when changing $\gamma =30$ to $\gamma =50$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ in at least one of its components for three vectors λ. For λ No. 3, we observe a deterioration of the first component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, and at the same time, an improvement of the second. For λ No. 5, we observe an improvement of the first component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, and at the same time, a deterioration of the second.
When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ.
For instance Bi9.1, we observe a similar phenomenon (an improvement, as well as a deterioration), although when changing $\gamma =10$ to $\gamma =30$, we observe an improvement in only one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for just two vectors λ. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for vectors λ Nos. 1–4 (at least one component improves). For λ No. 5, ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 $\gamma =10$ to $\gamma =30$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ in at least one of its components for three vectors λ, although only for one (λ No. 4) all components improve. When changing $\gamma =30$ to $\gamma =50$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ in at least one of its components for three vectors λ. We also observe a deterioration of the first component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for λ No. 2, and, at the same time, an improvement on its third component. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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_{\textit{opt}}}}}(\lambda ))$ are equal to the corresponding component of ${y^{\ast }}$. Only for λ No. 1, ${f_{3}}({x^{{P_{\textit{opt}}}}}(\lambda ))\lt {y_{3}^{\ast }}$, and when changing $\gamma =10$ to $\gamma =30$, as well as $\gamma =30$ to $\gamma =50$, there is an improvement only for the third component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$. 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 $\gamma =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 $\gamma =30$ and $\gamma =50$, the time to derive the corresponding upper shell increases significantly compared to $\gamma =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 ${T^{L}}+\text{Time}\hspace{0.1667em}{S_{U}}$.
5.5 Experiment 2 – Deriving Interval Representations with the Chute2 Algorithm
In this experiment, we check the behaviour of the Chute2 algorithm with $N=20$ and ${T^{S}}=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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Bi6.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Bi6.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Bi9.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | |||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Bi9.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma =10$ to $\gamma =30$, we observe an improvement in at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ, although only in three cases (λ Nos. 3–5) two components improve. Yet, when changing $\gamma =30$ to $\gamma =50$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, and, at the same time, an improvement on the first one. All components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ deteriorate for λ No. 4. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ (at least one component improves).
Table 22
Chute2, upper bounds for test problem Three6.1 and $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Three6.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Three6.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma \in \{10,30,50\}$.
| No. | $U({S_{U}},\lambda )$ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, ${\textit{GAP}_{{P_{sub}}}}\% $ for test problem Three9.1 and $\gamma \in \{10,30,50\}$.
| No. | ${\textit{GAP}_{{P_{sub}}}}\% $ | ||||||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =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, $|{S_{U}}|$, and Time$\hspace{0.1667em}{S_{U}}(s)$ for test problem Three9.1 and $\gamma \in \{10,30,50\}$.
| No. | $|{S_{U}}|$ | Time$\hspace{0.1667em}{S_{U}}(s)$ | ||||
| $\gamma =10$ | $\gamma =30$ | $\gamma =50$ | $\gamma =10$ | $\gamma =30$ | $\gamma =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 $\gamma =10$ to $\gamma =30$, we observe an improvement in at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ, and for λ Nos. 2–5, both components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improve. When changing $\gamma =30$ to $\gamma =50$, we observe an improvement in both components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ but the first one, where the first component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ deteriorates. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 $\gamma =10$, $\gamma =30$, and $\gamma =50$ are, respectively, 77.74, 85.03, and 83.28 seconds. These times are small fractions of the assumed time limit ${T^{L}}=1200$ seconds. For instance Bi9.1, average times over all vectors λ to derive the upper shell for $\gamma =10$, $\gamma =30$, and $\gamma =50$ are, respectively, 444.71, 535.92, and 618.98 seconds, and they are not small fractions of ${T^{L}}=1200$ seconds.
Let us check the results for tri-criteria instances. For instance Three6.1, when changing $\gamma =10$ to $\gamma =30$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, as well as a deterioration of the first one. When changing $\gamma =30$ to $\gamma =50$, we observe an improvement of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, as well as a deterioration of the first one. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in the ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for all vectors λ (at least one component improves).
For instance Three9.1, when changing $\gamma =10$ to $\gamma =30$, we observe an improvement in all components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for vectors λ Nos. 1–4, and for λ No. 5, ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ remains unchanged. When changing $\gamma =30$ to $\gamma =50$, we observe an improvement in all components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ for vector λ No. 1. For λ No. 2, we observe a deterioration of all components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, and for λ Nos. 3–4, we observe an improvement in two components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$, and a deterioration of one of its components. When changing $\gamma =10$ to $\gamma =50$, we observe an improvement in the ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 $\gamma =10$, $\gamma =30$, and $\gamma =50$ are, respectively, 425.48, 470.58, and 531.72 seconds, and they are not small fractions of ${T^{L}}=1200$ seconds.
We checked that for all instances, for no vector λ, the MIP solver derived the optimal solution to problem (2) within time limit ${T^{L}}+\text{Time}\hspace{0.1667em}{S_{U}}$.
5.6 Comparing 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$. We observe the following.
For instance Bi6.1, for all values of γ, all components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improve for all vectors λ.
For instance Bi9.1, for $\gamma =10$, all components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improve for four vectors λ, and one component improves for one vector λ. The same situation occurs for $\gamma =30$ and $\gamma =50$.
For instance Three6.1, for $\gamma =10$, at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improves for all vectors λ, and for two ones all components improve. The same situation occurs for $\gamma =30$ and $\gamma =50$.
For instance Three9.1, for all γ, at least one component of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improves for four vectors λ. All components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ improve for $\gamma =10$, $\gamma =30$, and $\gamma =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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$.
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 ${T^{S}}=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 ${T^{S}}$ 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 Time$\hspace{0.1667em}{S_{U}}(s)$ | AVG Time$\hspace{0.1667em}{S_{U}}(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.7 Discussion
For all test instances of the MOMIP problem, with time limits set, algorithm Chute2 determines tighter upper bounds measured with the help of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$. 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 ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ has been recorded for all vectors λ.
The deterioration of some of the components of ${G_{{P_{sub}}}}(R({S_{L}},{S_{U}},\lambda ))$ after increasing γ may be due to the fact that increasing the value of γ does not preserve the elements of the set ${S_{U}}$ obtained for smaller γ, but generates a new, denser set ${S_{U}}$, yet different in general. These new ${S_{U}}$ 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. ${T^{L}}=1200$ seconds, ${T^{S}}=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 ($\lambda ,\mu $) 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 $\Lambda :=\{{\lambda ^{1}},{\lambda ^{2}},{\lambda ^{3}}\}$ the algorithm derives upper shells ${S_{U}}({\lambda ^{1}})$, ${S_{U}}({\lambda ^{2}})$, and ${S_{U}}({\lambda ^{3}})$. ${S_{U}}:={\widehat{\oplus }_{i=1}^{3}}{S_{U}}({\lambda ^{i}})$ (where “$\widehat{\oplus }$” is an operator of adding two sets and removing dominating elements) is an upper shell, and ${S_{L}}:={\oplus _{i=1}^{3}}\{IN{C^{{\lambda ^{i}}}}\}$ (where “⊕” is an operator of adding two sets and removing dominated elements) is a lower shell. One can use ${S_{L}}$ and ${S_{U}}$ to calculate interval representations of implicit Pareto optimal outcomes designated by $\lambda \notin \Lambda $. For test problem Bi6.1 and $\gamma :=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 ${S_{L}}$ and upper shell ${S_{U}}$ elements in the objective space, respectively, ∙ – ${y^{\ast }}$.
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.
6 Limitations 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 ${\lambda ^{\prime }}$ 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 ${\lambda ^{\prime }}$ close to λ. However, this need not be the case, especially for vectors ${\lambda ^{\prime }}$, 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 ${\lambda ^{\prime }}$ considered in it (or, e.g. for ${\lambda ^{\prime }}$ 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 ${\lambda ^{\prime }}$ in the FindUpperShell algorithm. This time limit could be, e.g. a fraction of time ${T^{S}}$ adopted in Chute2. Since the number of probing vectors ${\lambda ^{\prime }}$ 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 ${\lambda ^{1}}=(0.055,0.945)$, $|{S_{U}}|=26$, Time$\hspace{0.1667em}{S_{U}}=91.44$ seconds, but for ${\lambda ^{5}}=(0.439,0.561)$, $|{S_{U}}|=7$, Time$\hspace{0.1667em}{S_{U}}=39.82$ seconds. To have the time of deriving an upper shell for some ${\lambda ^{\prime }}$ close to ${\lambda ^{1}}$ comparable to the time for ${\lambda ^{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 ${T^{L}}$, time limit ${T^{I}}$ (e.g. ${T^{L}}=1200$ seconds, ${T^{I}}=500$ seconds). Then the Chute algorithm would have time limit $\frac{{T^{I}}}{k}$ 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 ${T^{I}}$.
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 ${f_{l}}(x)\lt {y_{l}^{\ast }}$, $l=1,\dots ,k$. In this case, the upper bounds on components of Pareto optimal outcome designated by λ are not better than components of ${y_{l}^{\ast }}$. 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^{\ast }}$.
To find (sub)optimal values of multipliers ${\mu _{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).
7 Final Remarks
It has been shown how to algorithmically derive lower and upper shells to the MOMIP problem (for any $k\gt 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.
