∎ ∎

¹¹institutetext: W.T. Wu ²²institutetext: School of Sciences, Ningbo University of Technology 315211, China
²²email: weitianwu@nbut.edu.cn
🖂X.M. Yang ³³institutetext: National Center for Applied Mathematics of Chongqing 401331, China
School of Mathematical Sciences, Chongqing Normal University, Chongqing 401331, China
xmyang@cqnu.edu.cn

Obtaining properly Pareto optimal solutions of multiobjective optimization problems via the branch and bound method

Weitian Wu¹ Xinmin Yang²

(Received: date / Accepted: date)

Abstract

In multiobjective optimization, most branch and bound algorithms provide the decision maker with the whole Pareto front, and then decision maker could select a single solution finally. However, if the number of objectives is large, the number of candidate solutions may be also large, and it may be difficult for the decision maker to select the most interesting solution. As we argue in this paper, the most interesting solutions are the ones whose trade-offs are bounded. These solutions are usually known as the properly Pareto optimal solutions. We propose a branch-and-bound-based algorithm to provide the decision maker with so-called $\epsilon$ -properly Pareto optimal solutions. The discarding test of the algorithm adopts a dominance relation induced by a convex polyhedral cone instead of the common used Pareto dominance relation. In this way, the proposed algorithm excludes the subboxes which do not contain $\epsilon$ -properly Pareto optimal solution from further exploration. We establish the global convergence results of the proposed algorithm. Finally, the algorithm is applied to benchmark problems as well as to two real-world optimization problems.

Keywords:

Multiobjective optimization Global optimization Branch and bound algorithm Properly Pareto optimal solution

MSC:

90C2690C2990C30

1 Introduction

Many real-world optimization problems involve multiple objectives that require simultaneous consideration. Researchers refer to such problems as multiobjective optimization problems (MOPs). Since objectives are conflicting, it is impossible to find a single solution that is optimal for all objectives. Instead, a number of Pareto optimal solutions can be identified, characterized by the fact that improving one objective can only be achieved at the expense of worsening at least one other objective. Therefore, none of the Pareto optimal solutions can be said to be inferior when compared to others. The idea of solving an MOP can be understood as helping the decision maker in making a trade-off between multiple objectives and in finding a Pareto optimal solution that pleases him/her the most.

The a posteriori methods attempt to generate a well distributed set of representatives of the whole Pareto optimal solution set. Then the decision maker has to look at this potentially set of alternative solutions and make a choice. Most branch and bound algorithms for MOPs ref10 ; ref11 ; ref26 ; ref42 ; ref43 ; ref55 ; ref56 ; ref57 can be classified as a posteriori. However, their solution processes are considered resource-intensive and time-consuming, because the solutions with the required accuracy can only be obtained by continuously subdividing the variable space. Furthermore, if the number of objectives is large, not only the visualization of the high-dimensional Pareto front is not intuitive, but also the number of representative solutions may be huge. In this case, even if a global perspective of the problem is provided, it is very difficult to make a choice for decision maker.

To address the above issue, Wu and Yang ref44 proposed a reference point-based branch and bound algorithm which can provide the regions of interest that obey the decision maker’s preferences, instead of the whole Pareto front. The participation of preferences helps the algorithm to guide the search towards the regions of interest and avoiding computational resources being wasted on exploring undesirable regions. Thus, the algorithm not only requires less computation cost, but also effectively alleviates the decision maker’s selection pressure when solving many-objective optimization problems. However, the algorithm requires the assumption that the decision maker has informed preferences, i.e., the decision maker has sufficient a priori knowledge to specify preferences. This assumption does not hold in many situations ref21 , which makes the algorithm lose its advantages.

We note that some solutions are naturally preferred for the decision maker, even if no informed preference can be proposed. A classical type of the naturally preferred solutions is the properly Pareto optimal solutions. Kuhn and Tucker ref18 introduced the properly Pareto optimal solutions to avoid providing decision makers with some inappropriate Pareto optimal solutions whose trade-offs do not essentially differ from a weakly Pareto optimal solution. Therefore, from the view of trade-offs, the properly Pareto optimal solutions are more attractive to the decision maker than the improperly ones because their trade-offs are bounded. However, the methods available so far for finding the properly Pareto optimal solutions are almost based on the scalarization methodsref48 ; ref49 ; ref50 ; ref51 ; ref52 , so they cannot get rid of some limitations of scalarization methods. For instance, these methods require a set of predetermined parameters to control the distribution of the obtained solutions, however in some problems, choosing almost all parameters leads to unbounded problems ref53 . Furthermore, it is a frequent observation in some cases that even for convex problems, an evenly distributed set of parameters fails to produce an even distribution of solutions ref54 .

In this paper, we propose a branch-and-bound-based algorithm to approximate so-called $\epsilon$ -properly Pareto optimal solutions ref38 for nonconvex multiobjective optimization problems whose objective functions have the known Lipschitz constants. The discarding test of the proposed algorithm adopts a new dominance relation induced by a convex polyhedral ordering cone instead of the Pareto dominance relation. We prove the discarding test with the new dominance relation is able to exclude the subboxes which do not contain $\epsilon$ -properly Pareto optimal solutions, and further establish that the global convergence of the algorithm. Finally, the proposed algorithm is applied to two- and five-objective engineering constrained optimization problems as well as to benchmark problems. The proposed algorithm has the following capabilities:

•

Informed preferences of the decision maker are not required;
•

The algorithm approximates the $\epsilon$ -properly Pareto optimal solution set, instead of the Pareto optimal set;
•

The global convergence of the algorithm can be proved;
•

The algorithm is applicable to three or more objectives, and linear or nonlinear constraints.

The rest of this paper is organized as follows. In Section 2, we introduce the basic concepts and notations of multiobjective optimization. The new branch and bound algorithm and its theoretical analysis is described in Section 3. Section 4 is devoted to some numerical results.

2 Basics of multiobjective optimization

In this section we introduce the basic concepts which we need for the new algorithm. A multiobjective optimization problem can be written as follows:

\displaystyle\min\limits_{x\in\Omega}\quad F(x)=(f_{1}(x),\ldots,f_{m}(x))^{T}

(2.1)

with

\displaystyle\Omega=\{x\in\mathbb{R}^{n}:g_{j}(x)\geq 0,\;j=0,\ldots,p,\;\underline{x}_{k}\leq x_{k}\leq\overline{x}_{k},\;k=0,\ldots,n\},

where $f_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ ( $i=1,\ldots,m$ ) are Lipschitz continuous, and $g_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ ( $j=0,\ldots,p$ ) are continuous. If we allow $j=0$ , the set $\Omega$ is referred to as a box constraint. In this case, we call $\Omega$ a box with the midpoint $m(\Omega)=(\frac{\underline{x}_{1}+\overline{x}_{1}}{2},\ldots,\frac{\underline{x}_{n}+\overline{x}_{n}}{2})^{T}$ and the width $\omega(\Omega)=(\overline{x}_{1}-\underline{x}_{1},\ldots,\overline{x}_{n}-\underline{x}_{n})^{T}$ . The diameter of $\Omega$ is denoted by $\|\omega(\Omega)\|$ . For a feasible solution $x\in\Omega$ , the objective vector $F(x)\in\mathbb{R}^{m}$ is said to be the image of $x$ , while $x$ is called the preimage of $F(x)$ .

The concept of Pareto dominance relation between two solutions $x^{1},x^{2}\in\Omega$ can be defined as follows:

\displaystyle x^{1}\;dominates\;x^{2}~{}\Longleftrightarrow\;x^{1}\leq x^{2}\;\Longleftrightarrow\;F(x^{2})-F(x^{1})\in\mathbb{R}^{m}_{+}\backslash\{0\}.

We can define similar terms in the objective space. A nonempty set $N\subseteq\mathbb{R}^{m}$ is called a nondominated set if for any $z^{1},z^{2}\in N$ we have $z^{1}\nleq z^{2}$ and $z^{2}\nleq z^{1}$ .

A point $x^{*}\in\Omega$ is said to be a Pareto optimal solution for (MOP) if there does not exist any $x\in\Omega$ such that $F(x)\leq F(x^{*})$ . The set of all Pareto solutions is called the Pareto set and is denoted by $X^{*}$ . The image of Pareto set under the mapping $F$ is called the Pareto front.

The aim of an approximation algorithm is used to found an $\varepsilon$ -efficient solution of problem (2.1), which is defined next. Let $e$ denote the $m$ -dimensional all-ones vector $(1,\ldots,1)^{T}\in\mathbb{R}^{m}$ .

Definition 1.

ref47 Let $\varepsilon\geq 0$ be given. A point $\bar{x}\in\Omega$ is an $\varepsilon$ -efficient solution of problem (2.1) if there does not exist another $x\in\Omega$ with $F(x)\leq F(\bar{x})-\varepsilon e$ .

We use $d(a,b)=\|a-b\|$ to quantify the distance between two points $a$ and $b$ , where $\|\cdot\|$ denotes the Euclidean norm. The distance between the point $a$ and a nonempty finite set $B$ is defined as $d(a,B):=\min_{b\in B}\|a-b\|.$ Let $A$ be another non-empty finite set, we define the Hausdorff distance between $A$ and $B$ by

\displaystyle d_{H}(A,B):=\max\{d_{h}(A,B),d_{h}(B,A)\},

where $d_{h}(A,B)$ is the directed Hausdorff distance from $A$ to $B$ , defined by

\displaystyle d_{h}(A,B):=\max_{a\in A}\{\min_{b\in B}\|a-b\|\}.

3 Proposed branch and bound algorithm

Branch and bound algorithms ref10 ; ref11 ; ref26 ; ref42 ; ref43 ; ref55 ; ref56 ; ref57 have been used as the a posteriori methods to solve multiobjective optimization problems. By means of a tree search, a branch and bound algorithm systematically searches for an approximation of the entire Pareto set. The basic branch and bound algorithm for MOPs has been proposed by Fernández and Tóth ref11 (see Algorithm 1). The solution process consists of three components:

•

branching: subboxes are bisected perpendicularly to the direction of maximum width;
•

bounding: the lower and upper bounds for subboxes are calculated;
•

pruning: the subboxes that are provably suboptimal are excluded from exploration.

Input : problem (2.1), termination criterion;

Output :

\mathcal{B}_{k}

\mathcal{U}^{nds}

;

\mathcal{B}_{0}\leftarrow\Omega

\mathcal{U}^{nds}\leftarrow\emptyset

k=0

;

2 while termination criterion is not satisfied do

\mathcal{B}_{k+1}\leftarrow\emptyset

;

4 while $\mathcal{B}_{k}\neq\emptyset$ do

5 Select

B\in\mathcal{B}_{k}

and remove it from

\mathcal{B}_{k}

;

B_{1},B_{2}\longleftarrow

Bisect

B

perpendicularly to the direction of maximum width;

7 for $i=1,2$ do

8 Calculate the lower bound

l(B_{i})

and upper bound

u(B_{i})

for

B_{i}

;

9 if $B_{i}$ can not be discarded then

10 Update

\mathcal{U}^{nds}

u(B_{i})

and store

B_{i}

into

\mathcal{B}_{k+1}

;

11 end if

13 end for

15 end while

k\leftarrow k+1

17 end while

Algorithm 1 A basic branch and bound algorithm for MOPs

The source of the upper bounds is usually the image of the midpoints or the vertexes of subboxes. The approaches for the lower bounds proposed so far in the literature include the natural interval extension ref11 ; ref25 , the Lipschitz bound ref42 ; ref43 and the $\alpha$ BB method ref26 , and the resulting lower bound $l=(l_{1},\ldots,l_{m})^{T}$ for a subbox $B$ satisfies the following condition:

\displaystyle l\leq F(x),\quad x\in B.

(3.1)

Numerical experiments show that there is no big difference among the three bounding approaches. However, the latter two calculate the maximal error between the lower bounds and optimal values.

The pruning can be achieved by discarding tests. The discarding test limits the tree search and thus avoids exhaustive enumeration. A common type of discarding test is based on the Pareto dominance relation:

A subbox will be discarded if there exists a feasible objective vector such that the objective vector dominates the lower bound of the subbox.

It is easy to see that the subbox is removed because it does not contain any Pareto optimal solutions.

3.1 The new discarding test

As discussed earlier, we aim to approximate the properly Pareto optimal solutions. Here we consider the $\epsilon$ -properly Pareto optimal solution proposed by Wierzbicki ref38 :

Definition 2.

The solution $x^{*}\in\Omega$ is said to be the $\epsilon$ -properly Pareto optimal solution of problem (2.1), if

\displaystyle(F(x^{*})-\mathbb{R}^{m}_{\epsilon}\backslash\{0\})\cap F(\Omega)=\emptyset,

where $\mathbb{R}^{m}_{\epsilon}=\{y\in\mathbb{R}^{m}:\min_{i=1,\ldots,m}(1-\epsilon)y_{i}+\epsilon\sum_{i=1}^{m}y_{i}\geq 0\}$ , $0\leq\epsilon\leq 1$ .

Figure 1 depicts the $\epsilon$ -properly Pareto optimal solution of the bi-objective optimization problem. The $\epsilon$ -properly Pareto optimal solution can be obtained by intersecting the feasible region with a blunt cone $\mathbb{R}^{m}_{\epsilon}$ . Compared to the Pareto optimal solution, the $\epsilon$ -properly Pareto optimal solution uses a larger set $\mathbb{R}^{m}_{\epsilon}$ instead of $\mathbb{R}^{m}_{+}$ , so the $\epsilon$ -properly Pareto optimal solution set is contained in the Pareto optimal solution set. Furthermore, an interesting aspect of $\epsilon$ -properly Pareto optimal solutions is that the trade-offs are bounded by $\epsilon$ and $1/\epsilon$ ref23 ; ref59 .

Refer to caption — Figure 1: $\epsilon$ -properly Pareto optimal solution

For $0\leq\epsilon\leq 1$ , now we define a linear mapping $\mathcal{T}_{\epsilon}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ ,

\displaystyle\mathcal{T}_{\epsilon}(y):=\begin{bmatrix}1\quad&\quad\epsilon\quad&\quad\cdots\quad&\quad\epsilon\\ \epsilon\quad&\quad 1\quad&\quad\cdots\quad&\quad\epsilon\\ \vdots\quad&\quad\vdots\quad&\quad\ddots\quad&\quad\vdots\\ \epsilon\quad&\quad\epsilon\quad&\quad\cdots\quad&\quad 1\end{bmatrix}\cdot y.

Using this notation, we define a set

\displaystyle\mathcal{C}_{\epsilon}:=\{y\in\mathbb{R}^{m}:\mathcal{T}_{\epsilon}(y)\geqq 0\}.

By Definition 2.1.7 in ref58 , the set $\mathcal{C}_{\epsilon}$ is a convex polyhedral ordering cone, and thus the corresponding cone order $\leq_{\epsilon}$ on $\mathbb{R}^{m}$ by

\displaystyle y^{1}\leq_{\epsilon}y^{2}\Longleftrightarrow y^{2}-y^{1}\in\mathcal{C}_{\epsilon}\backslash\{0\}.

The following lemmas hold for $\leq_{\epsilon}$ :

Lemma 1.

ref13 For $y^{1},y^{2}\in\mathbb{R}^{m}$ , we have $y^{1}\leq_{\epsilon}y^{2}$ if and only if $\mathcal{T}_{\epsilon}(y^{1})\leq\mathcal{T}_{\epsilon}(y^{2})$ .

Proof.

we have

\displaystyle y^{1}\leq_{\epsilon}y^{2}\Leftrightarrow y^{2}-y^{1}\in\mathcal{C}_{\epsilon}\backslash\{0\}\Leftrightarrow\mathcal{T}_{\epsilon}(y^{2}-y^{1})\geq 0\Leftrightarrow\mathcal{T}_{\epsilon}(y^{2})\geq\mathcal{T}_{\epsilon}(y^{1}).

by the definition of $\leq_{\epsilon}$ and $\mathcal{C}_{\epsilon}$ and by the linearity of $\mathcal{T}_{\epsilon}$ .

Lemma 2.

The cone order $\leq_{\epsilon}$ is a strict partial order.

Proof.

By Definition 2.3.1 in ref58 , we would like to show that $\leq_{\epsilon}$ is irreflexive and transitive. First, it is easy to see that for all $y\in\mathbb{R}^{m}$ , $y\nleq_{\epsilon}y$ , i.e., $\leq_{\epsilon}$ is irreflexive. Next we will prove $\leq_{\epsilon}$ is transitive. For $y^{1},y^{2},y^{3}\in\mathbb{R}^{m}$ , assume $y^{1}\leq_{\epsilon}y^{2}$ and $y^{2}\leq_{\epsilon}y^{3}$ . By Lemma 1, we have $\mathcal{T}_{\epsilon}(y^{1})\leq\mathcal{T}_{\epsilon}(y^{2})$ and $\mathcal{T}_{\epsilon}(y^{2})\leq\mathcal{T}_{\epsilon}(y^{3})$ . Due to the transitivity of $\leq$ , we have $\mathcal{T}_{\epsilon}(y^{1})\leq\mathcal{T}_{\epsilon}(y^{3})$ , meaning that $y^{1}\leq_{\epsilon}y^{3}$ .

Lemma 3.

ref22 For $y^{1},y^{2}\in\mathbb{R}^{m}$ , if $y^{1}\leq y^{2}$ , we have $y^{1}\leq_{\epsilon}y^{2}$ .

Proof.

According to the Pareto dominance, we know $y^{2}-y^{1}\in\mathbb{R}^{m}_{+}\backslash\{0\}\subseteq\mathcal{C}_{\epsilon}\backslash\{0\}$ , following $\mathcal{T}_{\epsilon}(y^{2}-y^{1})\geq 0$ . By the linearity of $\mathcal{T}_{\epsilon}$ and Lemma 1, we have $y^{1}\leq_{\epsilon}y^{2}$ .

According to the strict partial order $\leq_{\epsilon}$ , the $\epsilon$ -dominance relation between $x^{1},x^{2}\in\Omega$ can be defined as follows:

\displaystyle x^{1}~{}\epsilon\hbox{-dominates}~{}x^{2}\Longleftrightarrow F(x^{1})\leq_{\epsilon}F(x^{2})\Longleftrightarrow\mathcal{T}_{\epsilon}(F(x^{1}))\leq\mathcal{T}_{\epsilon}(F(x^{2})).

A nonempty set $\mathcal{U}\subseteq\mathbb{R}^{m}$ is called non- $\epsilon$ -dominated set, if for any $y^{1}\in\mathcal{U}$ , there does not exist $y^{2}\in\mathcal{U}$ such that $y^{2}\leq_{\epsilon}y^{1}$ .

The following theorem discuss the relationship between $\epsilon$ -dominance and $\epsilon$ -properly Pareto optimal solutions.

Theorem 3.1.

For $\epsilon\in[0,1]$ , let $X^{*}_{\epsilon}$ be the $\epsilon$ -properly Pareto optimal solution set of problem (2.1). Then $x\in X^{*}_{\epsilon}$ if and only if there does not exist $x^{\prime}\in\Omega$ , such that $F(x^{\prime})\leq_{\epsilon}F(x)$ .

Proof.

( $\Rightarrow$ ) Assume that $x\in X^{*}_{\epsilon}$ . Suppose, to the contrary, that there exists $x^{\prime}\in\Omega$ , such that $F(x^{\prime})\leq_{\epsilon}F(x)$ . By Lemma 1, we have

\displaystyle F(x^{\prime})\leq_{\epsilon}F(x)\Leftrightarrow\mathcal{T}_{\epsilon}(F(x^{\prime}))\leq\mathcal{T}_{\epsilon}(F(x))\Leftrightarrow F(x)-F(x^{\prime})\in\mathcal{C}_{\epsilon}\backslash\{0\}.

Furthermore, it is easy to prove $\mathbb{R}^{m}_{\epsilon}=\mathcal{C}_{\epsilon}$ , thus we have $F(x)\in F(x^{\prime})+\mathbb{R}^{m}_{\epsilon}\backslash\{0\}$ , which is a contradiction to the fact that $x\in X^{*}_{\epsilon}$ .

( $\Leftarrow$ ) Assume that for $x\in\Omega$ , there does not exist $x^{\prime}\in\Omega$ , such that $F(x^{\prime})\leq_{\epsilon}F(x)$ and to the contrary $x\notin X^{*}_{\epsilon}$ . According to Definition 2, we know that there exists $\bar{x}\in\Omega$ , such that $F(x)\in F(\bar{x})+\mathcal{C}_{\epsilon}\backslash\{0\}$ . By Lemma 1 and the linearity of $\mathcal{T}_{\epsilon}$ , we obtain a contradiction $F(\bar{x})\leq_{\epsilon}F(x)$ .

According to Theorem 1, it is easy to see that the $\epsilon$ -properly Pareto optimal solution set can be found from a known Pareto optimal solution set by $\epsilon$ -dominance relation. On the other hand, the multiobjective branch and bound algorithms with the Pareto dominance-based discarding test can efficiently approximate the Pareto optimal solution set, thus the Pareto dominance in the discarding test can be generalized to the $\epsilon$ -dominance relation to ensure that $\epsilon$ -properly Pareto optimal solutions are identified. As a result, we can derive the following discarding test:

$\epsilon$ -Discarding test Let problem (2.1) be given, let $B$ be a subbox and $l(B)$ its lower bound. For $\epsilon\in[0,1]$ , if there exists a feasible objective vector $u\in F(\Omega)$ , such that $u\leq_{\epsilon}l(B)$ , then $B$ will be discarded.

Now we state the correctness of the proposed discarding test.

Theorem 3.2.

Let a subbox $B\in\Omega$ and its lower bound $l(B)\in\mathbb{R}^{m}$ be given, let $\mathcal{U}$ be a nondominated upper bound set of problem (2.1). For $\epsilon\in[0,1]$ , if there exists an upper bound $u\in\mathcal{U}$ such that $u\leq_{\epsilon}l(B)$ , then $B$ does not contain $\epsilon$ -properly Pareto optimal solution of problem (2.1).

Proof.

Let us assume that $x^{*}\in B$ is an $\epsilon$ -properly Pareto solution of problem (2.1). According to (3.1), we know that $l(B)\leq F(x^{*})$ . By Lemmas 2 and 3, we then have $u\leq_{\epsilon}l(B)\leq_{\epsilon}F(x^{*})$ . Thus from Theorem 1, we obtain a contradiction $x^{*}\notin X^{*}_{\epsilon}$ .

Algorithm 2 gives an implementation of the $\epsilon$ -discarding test, where the flag $D$ stands for decision to discard the subbox after the algorithm. Suppose that a non- $\epsilon$ -dominated upper bound set $\mathcal{U}$ is known, and if there exists an upper bound $u\in\mathcal{U}$ such that $u\leq_{\epsilon}l(B)$ , then the flag $D$ for $B$ is set to 1; otherwise, the flag $D$ is set to 0.

D\leftarrow 0

;

2 if there exists $u\in\mathcal{U}$ such that $u\leq_{\epsilon}l(B)$ then

D\leftarrow 1

;

5else

D\leftarrow 0

7 end if

return

D

Algorithm 2

\epsilon

-DiscardingTest(

B

l(B)

\mathcal{U}

)

3.2 The complete algorithm

Having the $\epsilon$ -discarding test, we can present the algorithm to find $\epsilon$ -properly Pareto solutions of problem (2.1). The whole algorithm is given in Algorithm 3. Another algorithmic design that needs to be mentioned is that we employ a parallel breadth first search strategy ref44 to search all subboxes simultaneously in each iteration. Therefore, in line 4 we simultaneously bisect all subboxes in the box collection $\mathcal{B}_{k-1}$ perpendicularly to the direction of maximum width to construct the current collection $\mathcal{B}_{k}$ . This branching produces subboxes that all have the same diameter, so in line 5 we only need to calculate the diameter of one subbox to obtain $\omega_{k}$ . After bisection, the feasibility test mentioned in ref11 will be used to exclude the subboxes that do not contain any feasible solutions from further exploration.

In the first for-loop, the lower bound $l(B)=(l_{1},\ldots,l_{m})^{T}$ of the subbox $B$ is calculated by

\displaystyle l_{i}=f_{i}(m(B))-\frac{L_{i}}{2}\|\omega(B)\|,\quad i=1,\ldots,m,

(3.2)

where $L_{i}$ is the Lipschitz constant of $f_{i}$ . Then, all lower bounds are stored in a lower bound set $\mathcal{L}$ . For each lower bound $l\in\mathcal{L}$ , we compare $l$ with other lower bounds stored in $\mathcal{L}$ by $\epsilon$ -dominance relation: if $l$ is $\epsilon$ -dominated by any other lower bounds, then $l$ will be removed from $\mathcal{L}$ ; otherwise, the lower bounds that are $\epsilon$ -dominated by $l$ will be removed from $\mathcal{L}$ . In this way, we can find a non- $\epsilon$ -dominated lower bound set $\mathcal{L}_{k}$ from $\mathcal{L}$ . The subboxes corresponding to the lower bounds stored in $\mathcal{L}_{k}$ constitutes $\bar{\mathcal{B}}$ .

In the second for-loop, an MOEA is used to find the upper bounds for all subboxes stored in $\bar{\mathcal{B}}$ . The upper bounds and the corresponding solutions (the preimage of the upper bounds) are stored in $\mathcal{U}$ and $\mathcal{X}$ , respectively. In order to reduce computational costs, we apply a mini MOEA which has a small initial population size and a few generations. If the problem also has the inequality constraints, the constrained handling approaches ref15 ; ref46 for MOEAs can be employed to ensure that feasible solutions and upper bounds are obtained. Thereafter, the non- $\epsilon$ -dominated upper bound set $\mathcal{U}_{k}$ and its corresponding solution set $\mathcal{X}_{k}$ will be found from $\mathcal{U}$ and $\mathcal{X}$ by means of $\epsilon$ -dominance relation.

In the third for-loop, the discarding flags for all subboxes are calculate by Algorithm 2 and are stored into a flag list $\mathcal{D}$ . Then, according to $\mathcal{D}$ , in line 22 the subboxes whose flags is equal to 1 will be removed from $\mathcal{B}_{k}$ . It is worth noting that, in order to speed up computation, we can compute the discarding flags for all the subboxes simultaneously, i.e., parallelize the third for-loop. The first and second for-loops can also be parallelized.

Input : problem (2.1),

\varepsilon>0

\delta>0

\epsilon\in[0,1]

;

Output :

\mathcal{B}_{k}

\mathcal{U}_{k}

\mathcal{X}_{k}

;

k\leftarrow 1

\mathcal{B}_{0}\leftarrow\Omega

\omega_{k-1}\leftarrow\|\omega(\Omega)\|

d\leftarrow 10^{6}

;

2 while $d>\varepsilon$ or $\omega_{k-1}>\delta$ do

\mathcal{L}\leftarrow\emptyset

\mathcal{U}\leftarrow\emptyset

\mathcal{X}\leftarrow\emptyset

;

4 Construct

\mathcal{B}_{k}

by bisecting all boxes in

\mathcal{B}_{k-1}

;

\omega_{k}\leftarrow\max\{\|\omega(B)\|:B\in\mathcal{B}_{k}\}

;

6 Update

\mathcal{B}_{k}

by the feasibility test suggested in ref11 ;

7 foreach $B\in\mathcal{B}_{k}$ do

8 Calculate for

B

its lower bound

l(B)

by equation (3.2);

\mathcal{L}\leftarrow\mathcal{L}\cup l(B)

;

11 end foreach

12 Find a non-

\epsilon

-dominated lower bound set

\mathcal{L}_{k}

from

\mathcal{L}

by the

\epsilon

-dominance;

13 Determine the box collection

\bar{\mathcal{B}}\subset\mathcal{B}_{k}

according to

\mathcal{L}_{k}

;

14 foreach $B\in\bar{\mathcal{B}}$ do

\bar{\mathcal{U}},\bar{\mathcal{X}}\leftarrow\textnormal{{MOEA}}(B)

;

\mathcal{U}\leftarrow\mathcal{U}\cup\bar{\mathcal{U}}

\mathcal{X}\leftarrow\mathcal{X}\cup\bar{\mathcal{X}}

;

18 end foreach

20 Find a non-

\epsilon

-dominated upper bound set

\mathcal{U}_{k}

and the corresponding solution set

\mathcal{X}_{k}

by the

\epsilon

-dominance;

21 foreach $B\in\mathcal{B}_{k}$ do

D(B)\leftarrow\textnormal{{$\epsilon$-DiscardingTest}}(B,l(B),\mathcal{U}_{k})

;

\mathcal{D}\leftarrow\mathcal{D}\cup D(B)

;

25 end foreach

26 Update

\mathcal{B}_{k}

according to the flag list

\mathcal{D}

;

d\leftarrow d_{h}(\mathcal{U}_{k},\mathcal{L}_{k})

k\leftarrow k+1

28 end while

Algorithm 3 Algorithm to find

\epsilon

-properly Pareto solutions

In order to ensure the trade-offs between disparately scaled objectives can be correctly expressed, we use an objective normalization technique:

\displaystyle\bar{f}_{i}=\frac{f_{i}-z^{*}_{i}}{z^{nad}_{i}-z^{*}_{i}},

where $z^{nad}=(z^{nad}_{1},\ldots,z^{nad}_{m})^{T}$ and $z^{*}=(z^{*}_{1},\ldots,z^{*}_{m})^{T}$ are the nadir and ideal points, respectively. There are two ways to obtain the nadir and ideal points: one is to solve each objective as a single-objective optimization problem to obtain the extreme points of the Pareto front. In this case, we should choose a global algorithm to solve the single-objective optimization problems; the second approach is to update $z^{nad}$ and $z^{*}$ during the solution process. The initial values of $z^{nad}$ and $z^{*}$ can be obtained from the natural interval extension of the objectives. The lower bound of the natural interval expansion is $z^{*}$ , and the upper bound is $z^{nad}$ . During the iterations, the minimum value of $f_{i}$ in the current lower bounds is used to update $z^{*}_{i}$ , while the maximum value of $f_{i}$ in the current upper bounds is used to update $\tilde{z}^{nad}$ . If a subbox can search for the maximum value of $f_{i}$ in the upper bounds or the minimum value of $f_{i}$ in the lower bounds, then it will be stored in the current box collection and not be removed by the $\epsilon$ -discarding test.

3.3 Convergence results

We start by showing the termination of the algorithm.

Theorem 3.3.

Let the predefined parameters $\varepsilon>0$ and $\delta>0$ be given, Algorithm 3 terminates after a finite number of iterations.

Proof.

Because we divide all boxes perpendicular to a side with maximal width, $\omega_{k}$ decreases among the sequence of box collections, i.e., $\omega_{k}>\omega_{k+1}$ for every $k$ and converges to 0. Therefore, for a given $\delta>0$ , there must exist a iteration count $\tilde{k}>0$ such that $\omega_{\tilde{k}}\leq\delta$ .

Assume we use (3.2) to calculate lower bounds. According to the way $\mathcal{U}_{k}$ is constructed and the $\epsilon$ -discarding test, for every $u\in\mathcal{U}_{k}$ , there exists a subbox $B\in\mathcal{B}_{k}$ , such that $F^{-1}(u)\in B$ and $l(B)\in\mathcal{L}_{k}$ . Then, based on the Lipschitz condition, we have

	$\displaystyle d(u,\mathcal{L}_{k})\leq d(u,l(B))=\\|u-F(m(B)+\frac{L}{2}\omega(B))\\|$	$\displaystyle\leq\\|u-F(m(B)\\|+\frac{1}{2}\omega_{k}\\|L\\|$
		$\displaystyle\leq\omega_{k}\\|L\\|$

where $L=(L_{1},\ldots,L_{m})^{T}$ consisting of the Lipschitz constants of objectives. Hence we have

\displaystyle d_{h}(\mathcal{U}_{k},\mathcal{L}_{k})=\omega_{k}\|L\|.

(3.3)

Due to the fact that $\omega_{k}$ converges to 0, it follows that for a given $\varepsilon>0$ , there must exist a iteration count $\bar{k}>0$ such that $d_{h}(\mathcal{U}_{k},\mathcal{L}_{k})\leq\varepsilon$ .

The next theorem states all properly Pareto optimal solutions of problem (2.1) are contained in the union of subboxes generated by the algorithm.

Theorem 3.4.

Proof.

Let us assume that there exists $k\in\mathbb{N}$ and an $\epsilon$ -properly Pareto solution $\bar{x}^{*}\in X_{\epsilon}^{*}$ such that $\bar{x}^{*}\notin\bigcup_{B\in\mathcal{B}_{k}}B$ , meaning that $\bar{x}^{*}$ is contained in a removed subbox $B$ in the pervious iteration. From the $\epsilon$ -discarding test, we know that $B$ will be discarded if and only if there exists a feasible objective vector $u$ such that $u\leq_{\epsilon}l(B)$ . Due to the lower bound (3.2) and Lemma 2, we then have $l(B)\leq_{\epsilon}F(\bar{x}^{*})$ . Thus we know that $u\leq_{\epsilon}F(\bar{x}^{*})$ , which contradicts the assumption that $\bar{x}^{*}$ is an $\epsilon$ -properly Pareto solution of problem (2.1).

Corollary 1.

Let $X_{\epsilon}^{*}$ be an $\epsilon$ -properly Pareto solution set of problem (2.1) and $\{\mathcal{B}_{k}\}_{k\in\mathbb{N}}$ a sequence of box collections generated by Algorithm 3. For each $x^{*}\in X_{\epsilon}^{*}$ and $k\in\mathbb{N}$ , these exists $B(x^{*},k)\in\mathcal{B}_{k}$ , such that $x^{*}\in B(x^{*},k)$ . Furthermore, we have $\lim\limits_{k\rightarrow\infty}d_{H}(X_{\epsilon}^{*},\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k))=0$ .

Proof.

The first conclusion is guaranteed by Theorem 4. Next we would like to prove the second conclusion.

On the one hand, from the first conclusion, we have

\displaystyle d(x^{*},\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k))=0,\quad x^{*}\in X^{*}_{\epsilon},~{}k\in\mathbb{N},

it follows $d_{h}(X^{*}_{\epsilon},\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k))=0$ . On the other hand, for each $x\in\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k)$ , there exists $B(\hat{x}^{*},k)\in\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k)$ , such that $x\in B(\hat{x}^{*},k)$ . We then have

\displaystyle 0\leq\lim\limits_{k\rightarrow\infty}d(x,X_{\epsilon}^{*})\leq\lim\limits_{k\rightarrow\infty}d(x,\hat{x}^{*})\leq\lim\limits_{k\rightarrow\infty}w_{k}=0,

meaning that $\lim\limits_{k\rightarrow\infty}d_{h}(\bigcup\limits_{x\in X_{\epsilon}^{*}}B(x,k),X_{\epsilon}^{*})=0$ . The second conclusion is proven.

In the next theorem, we prove that the solution set generated by Algorithm 3 is an $\varepsilon$ -efficient solution set of problem (2.1).

Theorem 3.5.

Assume that the midpoints of the subboxes are used to calculate the upper bounds in Algorithm 3. Let $\mathcal{X}$ be the solution set generated by Algorithm 3 and $\mathcal{L}_{k}$ the non- $\epsilon$ -dominated lower bound set. Then $\mathcal{X}$ is an $\varepsilon$ -efficient solution set of problem (2.1).

Proof.

Suppose $\tilde{x}\in\mathcal{X}$ is the midpoint of the subbox $B\in\mathcal{B}_{k}$ and $l\in\mathcal{L}_{k}$ is corresponding lower bound. According to (3.3), we know that

\displaystyle\varepsilon\geq\frac{1}{2}\omega_{k}\|L\|>\frac{1}{2}\omega_{k}L_{{\rm max}},

(3.4)

where $L_{max}=\max\{L_{i},i=1,\dots,m\}$ . Then we can obtain a lower bound $\tilde{l}=(\tilde{l}_{1},\ldots,\tilde{l}_{m})^{T}$ whose component can be calculated by

\displaystyle\tilde{l}_{i}=f(\tilde{x})_{i}-\frac{1}{2}\omega_{k}L_{{\rm max}},

and further, it is easy to see that $F(\tilde{x})-\varepsilon e<\tilde{l}\leqq l$ .

In the following we will prove $\tilde{x}$ is an $\varepsilon$ -efficient solution of problem (2.1) in two aspects. On the one hand, by (3.1), we have

\displaystyle F(\tilde{x})-\varepsilon e<\tilde{l}\leqq l\leq F(x),\quad x\in B.

Therefore, there does not exist $x\in B$ with $F(x)\leq F(\tilde{x})-\varepsilon e$ .

On the other hand, assume there exists another subbox $B^{\prime}\in\mathcal{B}_{k}\backslash B$ and a feasible point $x^{\prime}\in B^{\prime}$ such that $F(x^{\prime})\leq F(\tilde{x})-\varepsilon e$ . By (3.2), (3.4) and Lemma 3, we have

\displaystyle F(x^{\prime})-\frac{1}{2}\omega_{k}\|L\|\leq F(x^{\prime})\leq F(\tilde{x})-\varepsilon e<\tilde{l}\leqq l\Longleftrightarrow F(x^{\prime})-\frac{1}{2}\omega_{k}\|L\|\leq_{\epsilon}l,

which is a contradiction to the fact $\mathcal{L}_{k}$ is a non- $\epsilon$ -dominated lower bound set. Therefore, there does not exist a subbox $B^{\prime}\in\mathcal{B}_{k}\backslash B$ and a point $x^{\prime}\in B^{\prime}$ with $F(x^{\prime})\leq F(\tilde{x})-\varepsilon e$ . Now, we prove the conclusion.

To make it easier to obtain the conclusion in Theorem 3, we consider only the situation where the upper bounds are computed by using the midpoints, rather than mini MOEA. However, in practice, the tightness of the upper bounds obtained by mini MOEA tends to be better than the ones calculated by midpoints. Therefore, the accuracy of the solution set generated by Algorithm 3 will be higher than the accuracy of the one discussed in Theorem 5.

4 Experimental Results

Algorithm 3 is implemented in Python 3.8 with fundamental packages like numpy, scipy and multiprocessing, and performs on a computer with Intel(R) Core(TM) i7-10700 CPU and 32 Gbytes RAM on operation system WINDOWS 10 PROFESSIONAL. In all experiments, we set $\epsilon=0.75$ to search for the $\epsilon$ -properly Pareto solutions. For the MOEA employed in Algorithm 3, we use MOEA/D-DE ref20 with the population size 10 and 20 generations.

It is worth noting that we do not intend to compare Algorithm 3 with other algorithms. This is because the aim of most branch and bound algorithms is to search for Pareto optimal solution set rather than $\epsilon$ -properly Pareto solution set. Moreover, although some scalarization-based methods can obtain $\epsilon$ -properly Pareto solutions, not only they are not easily applied to some nonconvex problems, but also . Therefore, in the following we only demonstrate the effectiveness of Algorithm 3 on some test problems as well as engineering constrained optimization problems.

Example 1.

First, we consider three test problems.

•

MOP ref37

F(x)=\begin{pmatrix}0.5(\sqrt{1+(x_{1}+x_{2})^{2}}+\sqrt{1+(x_{1}-x_{2})^{2}}+x_{1}-x_{2})+e^{-(x_{1}-x_{2})^{2}}\\ 0.5(\sqrt{1+(x_{1}+x_{2})^{2}}+\sqrt{1+(x_{1}-x_{2})^{2}}-x_{1}+x_{2})+e^{-(x_{1}-x_{2})^{2}}\end{pmatrix},

where $x_{i}\in[-3,3],\;i=1,2$ . This problem has a disconnected Pareto front, which has two knees. We use $(\varepsilon,\delta)=(0.001,0.0001)$ .

•

DEB2DK ref3

F(x)=\begin{pmatrix}g(x)r(x_{1})\sin(0.5\pi x_{1})\\ g(x)r(x_{1})\cos(0.5\pi x_{1})\end{pmatrix},

where

	$\displaystyle g(x)=1+\frac{9}{n-1}\sum_{i=2}^{n}x_{i},$
	$\displaystyle r(x_{1})=5+10(x_{1}-0.5)^{2}+\frac{1}{K}\cos(2K\pi x_{1}),$
	$\displaystyle x_{i}\in[0,1],\;i=1,2.$

The parameter $K$ allows to control the number of knees in DEB2DK, and we set $K=4$ , $n=3$ and $(\varepsilon,\delta)=(0.0015,0.00015)$ .

•

DEB3DK ref3

F(x)=\begin{pmatrix}&g(x)r(x_{1},x_{2})\sin(0.5\pi x_{1})\sin(0.5\pi x_{2})\\ &g(x)r(x_{1},x_{2})\sin(0.5\pi x_{1})\cos(0.5\pi x_{2})\\ &g(x)r(x_{1},x_{2})\cos(0.5\pi x_{1})\end{pmatrix},

where

	$\displaystyle g(x)=1+\frac{9}{n-1}\sum_{i=3}^{n}x_{i},$
	$\displaystyle r(x_{1},x_{2})=(r_{1}(x_{1})+r_{2}(x_{2}))/2$
	$\displaystyle r_{i}(x_{i})=5+10(x_{i}-0.5)^{2}+\frac{1}{K}\cos(2K\pi x_{i}),$
	$\displaystyle x_{i}\in[0,1],\;i=1,2.$

In DEB3DK, we set $K=1$ , $n=3$ and $(\varepsilon,\delta)=(0.006,0.008)$ .

The experimental results for the three test problems are shown in Fig. 2. The Pareto fronts (blue dots) of three problems are found by Algorithm 3 with $\epsilon=0$ , and the red stars are the upper bounds obtained by Algorithm 3 with $\epsilon=0.75$ . An interesting observation can be made from the figure: the upper bounds obtained by Algorithm 3 usually lie in the “bulge” regions on Pareto fronts in Fig. 2. Based on this, we speculate that Algorithm 3 searches for the knees in the convex regions on the Pareto fronts. From the view of trade-offs, the convex knees on Pareto front are characterized by the fact that a small improvement in one objective will cause a large deterioration in the other objective ref3 . It is easy to see from Fig. 2 that $\epsilon$ -properly Pareto solutions found by Algorithm 3 also satisfy this feature. Existing knee-oriented approaches need to propose different indicators ref1 ; ref3 ; ref5 ; ref28 to identify knees, but these indicators are not actually mathematical definitions of knees. In contrast, $\epsilon$ -properly Pareto solutions not only satisfies the geometrical characterization of the knee, but also has a clear definition. Furthermore, several MOEAs ref22 ; ref29 also use $\epsilon$ -dominance relation or some similar dominance relation to find knees, but they do not discuss the relationship between the $\epsilon$ -dominance relation and $\epsilon$ -properly Pareto solutions, and thus there is no further proof of global convergence.

Example 2.

The welded beam design problem ref7 ; ref31 has four real-parameter variables $x=(x_{1},x_{2},x_{3},x_{4})$ and four non-linear constraints. Let us consider a manufacturing process in which minimization of the cost and minimization of the end deflection. The optimization problem is given as follows:

F(x)=\begin{pmatrix}1.10471x_{1}^{2}x_{3}+0.04811x_{2}x_{4}(14.0+x_{3})\\ 2.1592/(x_{2}x_{4}^{3})\end{pmatrix}

subject to the constraints

	$\displaystyle g_{1}(x)=13600-\tau(x)\geq 0,$
	$\displaystyle g_{2}(x)=30000-\sigma(x)\geq 0,$
	$\displaystyle g_{3}(x)=x_{2}-x_{1}\geq 0,$
	$\displaystyle g_{4}(x)=P_{c}(x)-6000\geq 0,$
	$\displaystyle 0.125\leq x_{1},x_{2}\leq 5,$
	$\displaystyle 0.1\leq x_{3},x_{4}\leq 10.$

where the stress and buckling terms are non-linear to design variables and are given as follows

	$\displaystyle\tau(x)=\sqrt{(\tau^{\prime})^{2}+(\tau^{\prime\prime})^{2}+(x_{3}\tau^{\prime}\tau^{\prime\prime})/\sqrt{0.25(x_{3}^{2}+(x_{1}+x_{4})^{2})}},$
	$\displaystyle\tau^{\prime}=\frac{6000}{\sqrt{2}x_{1}x_{3}},$
	$\displaystyle\tau^{\prime\prime}=\frac{6000(14+0.5x_{3})\sqrt{0.25(x_{3}^{2}+(x_{1}+x_{4})^{2})}}{1.414x_{1}x_{3}(x_{3}^{2}/12+0.25(x_{1}+x_{4})^{2})},$
	$\displaystyle\sigma(x)=\frac{504000}{x_{2}x_{4}^{2}},$
	$\displaystyle P_{c}(x)=64746.022(1-0.0282346x_{4})x_{4}x_{2}^{3}.$

Fig. LABEL:fig3 shows the results on the welded beam design problem. In this problem, we set $(\varepsilon,\delta)=(0.3,0.02)$ . The Pareto front are found by Algorithm 3 with $\epsilon=0$ . It can be seen that most of the objective vectors on the Pareto front have very high or very low trade-offs, meaning that these solutions do not differ from the weakly Pareto optimal solutions. Therefore, if the decision maker provides preferences without sufficient a priori information, he/she is likely to obtain undesired solutions with bad trade-offs. For example, we use the reference point-based branch and bound algorithm (RBB) mentioned in the literature ref44 to solve this problem. In RBB, we use the same $\varepsilon$ and $\delta$ , and assume that the decision maker provides three reference points as a priori information: (i) (4, 0.003), (ii) (20, 0.002), (iii) (32, 0.0007). The results of RBB are presented in Fig. 3b. It is not difficult to see that if the decision maker does not have a high demand for machining accuracy of the welded beam, they may not be satisfied with the solution corresponding to the second or third reference point, as he/she must pay high costs of fabrication to minimize the deflection. Therefore, the decision maker may adjust their preferences to the neighborhood of the first reference point. In contrast, Algorithm 3 directly provides the decision maker with $\epsilon$ -properly Pareto optimal solutions (red stars) without any a priori information. As can be seen in Fig. 3a, these solutions take into account both the cost and the deflection, and thus are more acceptable to the decision maker.

Example 3.

The water resource planning problem ref30 involves optimal planning for a storm drainage system in an urban area. The objectives to be minimized are drainage network cost, storage facility cost, treatment facility cost, expected flood damage cost and expected economic loss due to flood. The detailed description of the problem and constraints can be obtained from ref24 :

F(x)=\begin{pmatrix}106780.37(x_{2}+x_{3})+61704.67\\ 3000x_{1}\\ 30570*0.02289.0x_{2}/(0.06*2289.0)^{0.65}\\ 250.0*2289.0\exp(-39.75x_{2}+9.9x_{3}+2.74)\\ 25.0((1.39/(x_{1}x_{2}))+4940.0x_{3}-80.0)\end{pmatrix}

subject to the constraints

	$\displaystyle g_{1}(x)$	$\displaystyle=0.00139/(x_{1}x_{2})+4.94x3-0.08\leq 1,$
	$\displaystyle g_{2}(x)$	$\displaystyle=0.000306/(x_{1}x_{2})+1.082x_{3}-0.0986\leq 1,$
	$\displaystyle g_{3}(x)$	$\displaystyle=12.307/(x_{1}x_{2})+49408.24x_{3}+4051.02\leq 50000,$
	$\displaystyle g_{4}(x)$	$\displaystyle=2.098/(x_{1}x_{2})+8046.33x-696.71\leq 16000,$
	$\displaystyle g_{5}(x)$	$\displaystyle=2.138/(x_{1}x_{2})+7883.39x_{3}-705.04\leq 10000,$
	$\displaystyle g_{6}(x)$	$\displaystyle=0.417(x_{1}x_{2})+1721.26x3-136.54\leq 2000,$
	$\displaystyle g_{7}(x)$	$\displaystyle=0.164/(x_{1}x_{2})+631.13x3-54.48\leq 550,$
	$\displaystyle 0.01$	$\displaystyle\leq x_{1}\leq 0.45,~{}0.01\leq x_{2},x_{3}\leq 0.1.$

In this problem, we set $(\varepsilon,\delta)=(0.1,0.02)$ . The results of Algorithm 3 on the water resource planning problem are depicted in Fig. 4. The black lines are the value paths of the Pareto front obtained by Algorithm 3 with $\epsilon=0$ . Algorithm 3 with $\epsilon=0$ obtains a large number of candidates to represent the high-dimensional Pareto front, we only select some of them as representative solutions. The value paths of the images of $\epsilon$ -properly Pareto optimal solutions obtained by Algorithm 3 are plotted in red lines. It can be seen that the solutions provided by Algorithm 3 are only a small part of the entire Pareto front, so the decision maker does not confront a high level of decision pressure.

5 Conclusion

Many branch and bound algorithms for MOPs aim to approximate the whole Pareto optimal solution set. Their solution processes are considered resource-intensive and time-consuming. In particular, if the number of objectives is large, it is very difficult for the decision maker to select the most interesting solution from a large number of candidate solutions. Although the reference point-based branch and bound algorithm ref44 can obtain the regions of interest that match the decision maker’s preferences, it is not feasible to require the decision maker to explicitly specify a priori preferences in some situation. As a result, these algorithms may not be easy to use for decision makers.

In this paper, we have argued that, without the a priori preference, the properly Pareto optimal solutions are likely to be the most relevant to the decision maker. Consequently, we have proposed a new branch and bound algorithm to find so-called $\epsilon$ -properly Pareto optimal solutions. The basic idea was to replace the Pareto dominance relation in the discarding test with the $\epsilon$ -dominance relation that is induced by a convex polyhedral cone. In this way, the subboxes which do not contain the $\epsilon$ -properly Pareto optimal solution will be removed, resulting in a significant reduction in the number of candidate solutions. We prove that our algorithm is able to obtain a set of approximate solutions in a finite iteration. Numerical experiments confirm the effectiveness and application value of the proposed algorithm.

We are currently working on a refined version of the proposed algorithm, which allows to control the distribution of properly Pareto optimal solutions according to the skews of the Pareto front. Furthermore, it would be interesting to propose a branch and bound algorithm for vector optimization problems, and to test it on some real-world problems.

Acknowledgements.

This work is supported by the Major Program of National Natural Science Foundation of China (Nos. 11991020, 11991024), the General Program of National Natural Science Foundation of China (No. 11971084), the Team Project of Innovation Leading Talent in Chongqing (No. CQYC20210309536) and the NSFC-RGC (Hong Kong) Joint Research Program (No. 12261160365).

References

(1) Bhattacharjee, K. S., Singh, H. K., Ryan, M., & Ray, T. (2017). Bridging the gap: Many-objective optimization and informed decision-making. IEEE Transactions on Evolutionary Computation, 21(5), 813-820.
(2) Branke, J., Deb, K., Dierolf, H., & Osswald, M. (2004). Finding knees in multi-objective optimization. In Parallel Problem Solving from Nature-PPSN VIII: 8th International Conference, Birmingham, UK, September 18-22, 2004. Proceedings 8 (pp. 722-731). Springer Berlin Heidelberg.
(3) Cacchiani, V., & D¡¯Ambrosio, C. (2017). A branch-and-bound based heuristic algorithm for convex multi-objective MINLPs. European Journal of Operational Research, 260(3), 920-933.
(4) Das, I. (1999). On characterizing the ¡°knee¡± of the Pareto curve based on normal-boundary intersection. Structural optimization, 18, 107-115.
(5) De Santis, M., Eichfelder, G., Niebling, J., & Rocktäschel, S. (2020). Solving multiobjective mixed integer convex optimization problems. SIAM Journal on Optimization, 30(4), 3122-3145.
(6) Deb, K., & Sundar, J. (2006). Reference point based multi-objective optimization using evolutionary algorithms. In Proceedings of the 8th annual conference on Genetic and evolutionary computation (pp. 635-642).
(7) Eichfelder, G., Kirst, P., Meng, L., & Stein, O. (2021). A general branch-and-bound framework for continuous global multiobjective optimization. Journal of Global Optimization, 80, 195-227.
(8) Fernández, J., & Tóth, B. (2009). Obtaining the efficient set of nonlinear biobjective optimization problems via interval branch-and-bound methods. Computational Optimization and Applications, 42, 393-419.
(9) Fliege, J., Drummond, L. G., & Svaiter, B. F. (2009). Newton’s method for multiobjective optimization. SIAM Journal on Optimization, 20(2), 602-626.
(10) Ghaznavi, M., Akbari, F., & Khorram, E. (2021). Optimality conditions via a unified direction approach for (approximate) efficiency in multiobjective optimization. Optimization Methods and Software, 36(2-3), 627-652.
(11) Herzel, A., Helfrich, S., Ruzika, S., & Thielen, C. (2023). Approximating biobjective minimization problems using general ordering cones. Journal of Global Optimization, 86(2), 393-415.
(12) Hoseinpoor, N., & Ghaznavi, M. (2023). Finding properly efficient solutions of nonconvex multiobjective optimization problems with a minimum bound for trade-offs. International Journal of Nonlinear Analysis and Applications.
(13) Hozzar, B., Tohidi, G., & Daneshian, B. (2019). Two methods for determining properly effcient solutions with a minimum upper bound for trade-offs. Filomat, 33(6), 1551-1559.
(14) Jain, H., & Deb, K. (2013). An evolutionary many-objective optimization algorithm using reference-point based nondominated sorting approach, part II: Handling constraints and extending to an adaptive approach. IEEE Transactions on evolutionary computation, 18(4), 602-622.
(15) Jan, M. A., & Zhang, Q. (2010, September). MOEA/D for constrained multiobjective optimization: Some preliminary experimental results. In 2010 UK Workshop on computational intelligence (UKCI) (pp. 1-6). IEEE.
(16) Khaledian, K., & Soleimani-Damaneh, M. (2015). On efficient solutions with trade-offs bounded by given values. Numerical Functional Analysis and Optimization, 36(11), 1431-1447.
(17) Kuhn, H. W., & Tucker, A. W. (2013). Nonlinear programming. In Traces and emergence of nonlinear programming (pp. 247-258). Basel: Springer Basel.
(18) Kutateladze, S. S. (1979). Convex $\varepsilon$ -programming. Soviet Math. Dokl, 20(2).
(19) Li, H., & Zhang, Q. (2008). Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II. IEEE transactions on evolutionary computation, 13(2), 284-302.
(20) Li, K., Liao, M., Deb, K., Min, G., & Yao, X. (2020). Does preference always help? A holistic study on preference-based evolutionary multiobjective optimization using reference points. IEEE Transactions on Evolutionary Computation, 24(6), 1078-1096.
(21) Ikeda, K., Kita, H., & Kobayashi, S. (2001). Failure of Pareto-based MOEAs: Does non-dominated really mean near to optimal?. In Proceedings of the 2001 congress on evolutionary computation (IEEE Cat. No. 01TH8546) (Vol. 2, pp. 957-962). IEEE.
(22) Miettinen, K. (1999). Nonlinear multiobjective optimization. Springer Science & Business Media.
(23) Morovati, V., & Pourkarimi, L. (2019). Extension of Zoutendijk method for solving constrained multiobjective optimization problems. European Journal of Operational Research, 273(1), 44-57.
(24) Musselman, K., & Talavage, J. (1980). A tradeoff cut approach to multiple objective optimization. Operations Research, 28(6), 1424-1435.
(25) Neumaier, A. (1990). Interval methods for systems of equations (No. 37). Cambridge university press.
(26) Niebling, J., & Eichfelder, G. (2019). A branch–and–bound-based algorithm for nonconvex multiobjective optimization. SIAM Journal on Optimization, 29(1), 794-821.
(27) Przybylski, A., & Gandibleux, X. (2017). Multi-objective branch and bound. European Journal of Operational Research, 260(3), 856-872.
(28) Rachmawati, L., & Srinivasan, D. (2009). Multiobjective evolutionary algorithm with controllable focus on the knees of the Pareto front. IEEE Transactions on Evolutionary Computation, 13(4), 810-824.
(29) Ramirez-Atencia, C., Mostaghim, S., & Camacho, D. (2017). A knee point based evolutionary multi-objective optimization for mission planning problems. In Proceedings of the genetic and evolutionary computation conference (pp. 1216-1223).
(30) Ray, T., Tai, K., & Seow, K. C. (2001). Multiobjective design optimization by an evolutionary algorithm. Engineering Optimization, 33(4), 399-424.
(31) Ravindran, A., Reklaitis, G. V., & Ragsdell, K. M. (2006). Engineering optimization: methods and applications. John Wiley & Sons.
(32) Sawaragi, Y., Nakayama, H., & Tanino, T. (1985). Theory of multiobjective optimization. Elsevier.
(33) Schütze, O., Laumanns, M., & Coello, C. A. C. (2008). Approximating the knee of an MOP with stochastic search algorithms. In Parallel Problem Solving from Nature¨CPPSN X: 10th International Conference, Dortmund, Germany, September 13-17, 2008. Proceedings 10 (pp. 795-804). Springer Berlin Heidelberg.
(34) Tang, L., & Yang, X. (2021). A modified direction approach for proper efficiency of multiobjective optimization. Optimization Methods and Software, 36(2-3), 653-668.
(35) Wierzbicki, A. P. (1980). The use of reference objectives in multiobjective optimization. In Multiple Criteria Decision Making Theory and Application: Proceedings of the Third Conference Hagen/Königswinter, West Germany, August 20¨C24, 1979 (pp. 468-486). Springer Berlin Heidelberg.
(36) Wierzbicki, A. P. (1986). A methodological approach to comparing parametric characterizations of efficient solutions. In Large-Scale Modelling and Interactive Decision Analysis: Proceedings of a Workshop sponsored by IIASA (International Institute for Applied Systems Analysis) and the Institute for Informatics of the Academy of Sciences of the GDR Held at the Wartburg Castle, Eisenach, GDR, November 18¨C21, 1985 (pp. 27-45). Springer Berlin Heidelberg.
(37) Wu, W. T., & Yang, X. M. (2023). Reference-point-based branch and bound algorithm for multiobjective optimization. Journal of Global Optimization, 1-19.
(38) Žilinskas, A., & Gimbutienė, G. (2016). On one-step worst-case optimal trisection in univariate bi-objective Lipschitz optimization. Communications in Nonlinear Science and Numerical Simulation, 35, 123-136.
(39) Žilinskas, A., & Žilinskas, J. (2015). Adaptation of a one-step worst-case optimal univariate algorithm of bi-objective Lipschitz optimization to multidimensional problems. Communications in Nonlinear Science and Numerical Simulation, 21(1-3), 89-98.