The Pareto Tracer for General Inequality Constrained Multi-Objective Optimization Problems

Abstract

Problems where several incommensurable objectives have to be optimized concurrently arise in many engineering and financial applications. Continuation methods for the treatment of such multi-objective optimization methods (MOPs) are very efficient if all objectives are continuous since in that case one can expect that the solution set forms at least locally a manifold. Recently, the Pareto Tracer (PT) has been proposed, which is such a multi-objective continuation method. While the method works reliably for MOPs with box and equality constraints, no strategy has been proposed yet to adequately treat general inequalities, which we address in this work. We formulate the extension of the PT and present numerical results on some selected benchmark problems. The results indicate that the new method can indeed handle general MOPs, which greatly enhances its applicability.

1. Introduction

In many real-world applications, the problem occurs that several conflicting and incommensurable objectives have to be optimized concurrently. As general example, in the design of basically any product, both cost (to be minimized) and quality (to be maximized) are relevant objectives, among others. Problems of that kind are termed multi-objective optimization problems (MOPs). In the case all of the objectives are continuous and in conflict with each other, it is known that there is not one single solution to be expected (as it is the case for scalar optimization problems, i.e., problems where one objective is considered) but an entire set of solutions. More precisely, one can expect that the solution set—the Pareto set, and, respectively its image, the Pareto front—forms at least locally an object of dimension $k-1$, where k is the number of objectives involved in the problem. Due to this, “curse of dimensionality” problems with more than, e.g., four objectives are also called many objective optimization problems (MaOPs).

In the literature, many different methods for the numerical treatment of MOPs and MaOPs can be found (see also the discussion in the next section). One class of such methods is given by specialized continuation methods that take advantage of the fact that the solution set forms—at least locally and under certain mild assumption on the model as discussed in [1]—a manifold. Continuation methods start with one (approximate) solution of the problem and perform a movement along the Pareto set/front of the given M(a)OP via considering the underdetermined system of equations that is developed out of the Karush–Kuhn–Tucker (KKT) equations of the problem. By construction, continuation methods are of local nature. That is, if the Pareto set consists of different connected components, such methods will have to be fed with several starting points in order to obtain approximations of the entire solution set. On the other hand, continuation methods are probably most effective locally (i.e., within each connected component). Thus far, several multi-objective continuation methods have been proposed. Most of these continuation methods, however, are designed for or restricted to the treatment of bi-objective problems (i.e., MOPs with two objectives). The method of Hillermeier [1] and the Pareto Tracer (PT [2]) have been proposed for general number k of objectives. The method of Hillermeier is applicable to unconstrained and equality constrained MOPs, and the PT in addition to box constrained problems. Thus far, no extensions for these two methods are known for the treatment of general inequalities, which represents a significant shortcoming since such constraints naturally arise in many applications (e.g., [3,4]). In this paper, we extend the PT for the treatment of general inequality constraints. To this end, we utilize and adapt elements from active set methods to decide which of the inequalities have to be treated as equalities at each candidate solution. We demonstrate the strength of the novel algorithm on several benchmark test functions and present comparisons to some other numerical multi-objective solvers. The results indicate that the new method can indeed reliably handle MOPs with general constraints.

The remainder of this paper is organized as follows. In Section 2, we shortly present the required background for the understanding of this work. In Section 3, we adapt the Pareto Tracer for the treatment of general (equality and inequality) constraints. In Section 4, we present some results of the PT as well as some other multi-objective numerical methods on selected benchmark problems. Finally, we draw our conclusions in Section 5 and give possible paths for future research.

2. Background and Related Work

In this section, we briefly state the main concepts and notations that are used for the understanding of this work (for details, we refer to, e.g., [5,6]).

We consider here continuous multi-objective optimization problem (MOPs) that can be defined mathematically as

$$\begin{array}{cc}\underset{x}{min}& F\left(x\right),\hfill \\ s.t.& h_i\left(x\right)=0,i=1,\dots ,p,\hfill \\ & g_i\left(x\right)\le 0,i=1,\dots ,m,\hfill \end{array}$$

(1)

where $F:Q\subset {\mathbb{R}}^n\to {\mathbb{R}}^k$, $F\left(x\right)={(f_1\left(x\right),\dots ,f_k\left(x\right))}^T$ is the map of the k individual objectives $f_i:Q\subset {\mathbb{R}}^n\to \mathbb{R}$. We assume that all objectives and constraint functions are twice continuously differentiable. The domain Q of the functions is defined by the equality and inequality constraints of (1):

$$Q:=\{x\in {\mathbb{R}}^n:h_i\left(x\right)=0,i=1,\dots ,p\mathrm{and}{g}_i\left(x\right)\le 0,i=1,\dots ,m\}.$$

(2)

If a point $x\in {\mathbb{R}}^n$ satisfies all constraints of (1), i.e., if $x\in Q$, we call this point feasible. Points $x\notin Q$ are called infeasible. If $k=2$ objectives are considered, the problem is also termed a bi-objective optimization problem (BOP).

We say that a point $x\in Q$ dominates a point $y\in Q$ (in short: $x\prec y$) if $f_i\left(x\right)\le f_i\left(y\right)$ for all $i=1,\dots ,k$, and there exists an index j such that $f_j\left(x\right)<f_j\left(y\right)$. A point $x^{\ast }$ is called Pareto optimal or simply optimal if there does not exist a vector $y\in Q$ that dominates $x^{\ast }$. A point $x^{\ast }\in Q$ is called locally optimal if there does not exist a vector $y\in Q\cap N\left(x^{\ast }\right)$ that dominates $x^{\ast }$, where $N\left(x^{\ast }\right)$ is a neighborhood of $x^{\ast }$. The set $P_Q$ of all Pareto optimal solutions is called the Pareto set, and its image $F\left(P_Q\right)$ the Pareto front. In [1], it has been shown that one can expect that both Pareto set and front typically form ($k-1$)-dimensional objects under certain (mild) conditions on the problem.

If all objectives and constraint functions are differentiable, local optimal solutions can be characterized by the Karush–Kuhn–Tucker (KKT) equations [7,8]:

Theorem 1.

Suppose that $x^{\ast }$ is locally optimal with respect to (1). Then, there exist Lagrange multipliers $\alpha \in {\mathbb{R}}^k$, $\lambda \in {\mathbb{R}}^p$ and $\gamma \in {\mathbb{R}}^m$ such that the following conditions are satisfied

$$\begin{array}{cc}\hfill \sum _{i=1}^k{\alpha }_i\nabla f_i\left(x^{\ast }\right)+\sum _{i=1}^p{\lambda }_i\nabla h_i\left(x^{\ast }\right)+\sum _{i=1}^m{\gamma }_i{g}_i\left(x^{\ast }\right)=0& \hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill h_i\left(x^{\ast }\right)=0,& i=1\dots p,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill g_i\left(x^{\ast }\right)\le 0,& i=1\dots m,\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\alpha }_i\ge 0,& i=1\dots k,\hfill \end{array}$$

(3d)

$$\begin{array}{cc}\hfill \sum _{i=1}^k{\alpha }_i=1,& \hfill \end{array}$$

(3e)

$$\begin{array}{cc}\hfill {\gamma }_i\ge 0,& i=1\dots m,\hfill \end{array}$$

(3f)

$$\begin{array}{cc}\hfill {\gamma }_i{g}_i\left(x^{\ast }\right)=0,& i=1\dots m.\hfill \end{array}$$

(3g)

Multi-objective optimization is an active field of research, and thus far many numerical methods have been proposed for the treatment of such problems. There exist for instance many methods that are designed to compute single solutions such as the weighted sum method [9], the $ϵ$-constraint method [5,10], the weighted metric and weighted Tchebycheff method [5,11,12], as well as reference point problems [13,14,15]. All of these methods transform the given MOP into a scalar optimization problem (SOP) that can to a certain extent to include users’ preferences. These methods can either be used as standalone algorithm (i.e., for the computation of single solutions) or be used to obtain a finite size approximation of the entire Pareto set/front of the given MOP via utilizing a clever sequence of these SOPs [5,16,17,18,19].

Further, there exist set oriented methods such as cell mapping techniques [20,21,22,23]), subdivision techniques [24,25,26,27], and multi-objective evolutionary algorithms (MOEAs, e.g., [3,28,29,30,31,32,33,34]). All of these methods manipulate an entire set of candidate solutions in each iteration and hence yield a finite size approximation of the solution set in one run of the algorithm. Hybridizations of such techniques with mathematical programming techniques can be found in [31,35,36,37,38,39,40,41].

Finally, a third class of numerical solvers for MOPs is given by specialized continuation methods that take advantage of the fact that the Pareto set/front of a given problem forms at least locally a manifold of a certain dimension. Methods of this kind start with a given (approximate) solution and perform a movement along the Pareto set/front of the problem. The first such method is proposed in [1], which can be applied to unconstrained and equality constrained MOPs of any number k of objectives, while no strategies are reported on how to treat inequalities. ParCont [42,43] is a rigorous predictor–corrector method that is based on interval analysis and parallelotope domains. The method can deal with equality and inequality constraints, but it is restricted to bi-objective problems. This restriction also holds for the method presented in [44], which has been designed to provide an equispaced approximation of the Pareto front. The Zigzag method [45,46,47] obtains Pareto front approximations via alternating optimizing one of the objectives. This approach is also limited to the treatment of bi-objective problems.

In [48], a continuation method is presented that is applicable to box-constrained BOPs. In [49], a variant of the method of Hillermeier is presented that is designed for the treatment of high-dimensional problems.

Recently, the Pareto Tracer (PT) was proposed by Martin and Schütze [2]. Similar to the method of Hillermeier, PT addresses the underdetermined nonlinear system of equations that is induced by the KKT equations. However, unlike the method of Hillermeier, the PT aims to separate the decision variables from the associated weight (or Lagrange) vectors whenever possible, leading to significant changes. The latter is due to the fact that the nonlinearity of the equation system can be significantly higher in the compound space compared to the corresponding system that is only defined in decision variable space. As a by-product, the chosen approach allows to compute the tangent space of both Pareto set and front at every given regular point x. In [50], elements of the PT are used to treat many objective optimization problems (i.e., MOPs with more than, e.g., five objectives). Thus far, PT is only applicable to box and equality constrained problems which limits its application. In the following, we propose and discuss an extension of this method to adequately treat general MOPs, i.e., MOPs that in particular contain general inequalities.

3. Adapting the Pareto Tracer for General Inequality Constrained MOPs

In this section, we adapt the PT so that is can handle general inequality constraints. The core is the predictor–corrector step that generates from a given candidate solution $x_i$ the following candidate $x_{i+1}$ that satisfies the KKT conditions, and so that $F\left(x_{i+1}\right)-F\left(x_i\right)$ defines a pre-described movement in objective space along the set of KKT points.

Assume we are given a MOP of form (1) and a feasible point $x_0$ that satisfies the KKT conditions (3), where ${\alpha }_i>0$, $i=1,\dots ,k$. Let $ϵ>0$ and define by

$$I_p\left(ϵ\right):=\{j\in \{1,\dots ,m\}:g_j\left(x_0\right)\ge -ϵ\}$$

(4)

the set of indices corresponding to the nearly active inequalities at $x_0$. If $I_p\left(ϵ\right)=\{j_1,\dots ,j_s\}$, $s\le m$, define

$$G_{ϵ}:=\left(\begin{array}{c}\hfill \nabla g_{j_1}{\left(x_0\right)}^T\\ ⋮\\ \hfill \nabla g_{j_s}{\left(x_0\right)}^T\end{array}\right)\in {\mathbb{R}}^{s\times n}.$$

(5)

Further, let

$$\begin{array}{cc}\hfill J& :=\left(\begin{array}{c}\hfill \nabla f_1{\left(x\right)}^T\\ ⋮\\ \hfill \nabla f_k{\left(x\right)}^T\end{array}\right)\in {\mathbb{R}}^{k\times n}\hfill \\ \hfill H& :=\left(\begin{array}{c}\hfill \nabla h_1{\left(x_0\right)}^T\\ ⋮\\ \hfill \nabla h_p{\left(x_0\right)}^T\end{array}\right)\in {\mathbb{R}}^{p\times n},\hfill \end{array}$$

(6)

and $\alpha \in {\mathbb{R}}^n$, $\lambda \in {\mathbb{R}}^p$, and $\gamma \in {\mathbb{R}}^s$ be the solution of

$$\underset{\overline{\alpha },\overline{\lambda },\overline{\gamma }}{min}\left\{{∥J^T\overline{\alpha }+H^T\overline{\lambda }+G_{ϵ}^T\overline{\gamma }∥}_2^2:{\overline{\alpha }}_i\ge 0,i=1,\dots ,k,\sum _{i=1}^k{\overline{\alpha }}_i=1\right\}.$$

(7)

Note that (7) yields the Lagrange multipliers at $x_0$ for $ϵ=0$ if $x_0$ is a KKT point and if all active inequalities are regarded as equalities. Using $\alpha$, $\lambda$ and $\gamma$, define the matrix

$$W_{\alpha ,\beta ,\gamma }:=\sum _{i=1}^k{\alpha }_i{\nabla }^2{f}_i\left(x\right)+\sum _{i=1}^p{\lambda }_i{\nabla }^2{h}_i\left(x\right)+\sum _{i=1}^s{\gamma }_i{\nabla }^2{g}_{j_i}\left(x\right)\in {\mathbb{R}}^{n\times n}.$$

(8)

To compute a predictor direction ${\nu }_{\mu }\in {\mathbb{R}}^n$, we solve the the system

$$\left(\begin{array}{ccc}{W}_{\alpha ,\lambda ,\gamma }& H^T& G_{ϵ}^T\\ H& 0& 0\\ G_{ϵ}& 0& 0\end{array}\right)\left(\begin{array}{c}{\nu }_{\mu }\\ \zeta \\ \sigma \end{array}\right)=\left(\begin{array}{c}-J^T\mu \\ 0\\ 0\end{array}\right).$$

(9)

Note that system (9) depends on $\mu \in {\mathbb{R}}^k$. Before we specify this vector, we first simplify (9). Denote by

$$A:=\left(\begin{array}{c}H\\ G_{ϵ}\end{array}\right)\in {\mathbb{R}}^{(p+s)\times n},\xi :=\left(\begin{array}{c}\zeta \\ \sigma \end{array}\right)\in {\mathbb{R}}^{p+s},$$

(10)

then (9) is equivalent to

$$\left(\begin{array}{cc}{W}_{\alpha ,\lambda ,\gamma }& A^T\\ A& 0\end{array}\right)\left(\begin{array}{c}{\nu }_{\mu }\\ \xi \end{array}\right)=\left(\begin{array}{c}-J^T\mu \\ 0\end{array}\right).$$

(11)

Let $d\in {\mathbb{R}}^k$. It is straightforward to show that for a vector ${\nu }_{{\mu }_d}$ that solves (11), where ${\mu }_d\in {\mathbb{R}}^k$ is chosen such that

$$\left(\begin{array}{c}-J{W}_{\alpha ,\lambda ,\gamma }^{-1}{J}^T\\ 1\dots 1\end{array}\right){\mu }_d=\left(\begin{array}{c}d\\ 0\end{array}\right),$$

(12)

it holds

$$J{\nu }_{{\mu }_d}=d.$$

(13)

That is, (infinitesimal) small steps from $x_0$ into direction ${\nu }_{{\mu }_d}$ (in decision variable space) will lead to a movement from $F\left(x_0\right)$ into direction d (in objective space). It remains to select a suitable choice for d. Since $\alpha$ is orthogonal to the linearized Pareto front at $F\left(x_0\right)$ [1], a suggesting choice is hence by (13) to take d orthogonal to $\alpha$. For this, let

$$\alpha =QR=(q_1,q_2,....,q_k)R,$$

(14)

where $Q\in {\mathbb{R}}^{k\times k}$ is orthogonal and $R\in {\mathbb{R}}^{k\times 1}$, be a $QR$-factorization of $\alpha$. Then, any vector

$$d\in span\{q_2,\dots ,q_k\}$$

(15)

can be chosen so that a movement in direction ${\nu }_{{\mu }_d}$ (in decision variable space) leads to a movement from $F\left(x_0\right)$ along the Pareto front. Note that the second equation in (12) reads as ${\sum }_{i=1}^k{\mu }_i=0$. Hence, for the special case of a bi-objective optimization problem (i.e., $k=2$), there are—after normalization—only two choices for $\mu$:

$${\mu }^{\left(1\right)}=\left(\begin{array}{c}-1\\ 1\end{array}\right),\mathrm{and}{\mu }^{\left(2\right)}=\left(\begin{array}{c}1\\ -1\end{array}\right).$$

(16)

Analog to Martin and Schütze [2], one can show that ${\mu }^{\left(1\right)}$ corresponds to a “right down” movement along the Pareto front while ${\mu }^{\left(2\right)}$ corresponds to a "left up" movement along the Pareto front.

After selecting the predictor direction ${\nu }_{\mu }$, the question is how far to step in this direction. Here, we follow the suggestion made by Hillermeier [1] and use the step size

$$t=\frac{\tau }{\parallel J{\nu }_{\mu }{\parallel }_2}$$

(17)

for a (small) value $\tau >0$ so that

$$\parallel F(x_0+t{\nu }_{\mu })-F\left(x_0\right){\parallel }_2\approx \tau .$$

(18)

For the computations presented below, we make the following modifications: instead of $W_{\alpha ,\beta ,\gamma }$, we use the matrix

$$W_{\alpha }:=\sum _{i=1}^k{\alpha }_i{\nabla }^2{f}_i\left(x\right)\in {\mathbb{R}}^{n\times n}.$$

(19)

More precisely, for the computation of ${\nu }_{\mu }$, we use the system

$$\left(\begin{array}{cc}{W}_{\alpha }& A^T\\ A& 0\end{array}\right)\left(\begin{array}{c}{\nu }_{\mu }\\ \xi \end{array}\right)=\left(\begin{array}{c}-J^T\mu \\ 0\end{array}\right)$$

(20)

and to obtain ${\mu }_d$ we solve

$$\left(\begin{array}{c}-J{W}_{\alpha }^{-1}{J}^T\\ 1\dots 1\end{array}\right){\mu }_d=\left(\begin{array}{c}d\\ 0\end{array}\right).$$

(21)

We have observed similar performance for both approaches, while the usage of $W_{\alpha }$ compared to $W_{\alpha ,\beta ,\gamma }$ comes with the advantage that no Hessians for any of the constraint functions have to be computed.

Given a predictor point

$${\tilde{x}}_1:=x_0+t{\nu }_{\mu }$$

(22)

the task of the upcoming corrector step is to find a KKT point $x_1$ that is ideally near to ${\tilde{x}}_1$. For this, we suggest to apply the multi-objective Newton method proposed in [51]. In particular, we first compute the solution $({\tilde{\nu }}_1,\tilde{\delta })$ of the following problem

$$\begin{array}{cc}\hfill \underset{(\nu ,\delta )\in {\mathbb{R}}^n\times \mathbb{R}}{min}\delta & \\ \hfill \mathrm{s}.\mathrm{t}.& \nabla f_i{\left({\tilde{x}}_1\right)}^T\nu +\frac{1}{2}{\nu }^T{\nabla }^2{f}_i\left({\tilde{x}}_1\right)\nu \le \delta ,i=1,...,k,\hfill \\ \hfill & h_i\left({\tilde{x}}_1\right)+\nabla h_i{\left({\tilde{x}}_1\right)}^T\nu =0,i=1,...,p.\hfill \end{array}$$

(23)

${\tilde{\nu }}_1$ is indeed the Newton direction for equality constrained MOPs as suggested in [2]. To adequately treat the involved inequalities, however, we propose to use the solution of the following problem:

$$\begin{array}{cc}\hfill \underset{(\nu ,\delta )\in {\mathbb{R}}^n\times \mathbb{R}}{min}\delta & \\ \hfill \mathrm{s}.\mathrm{t}.& \nabla f_i{\left({\tilde{x}}_1\right)}^T\nu +\frac{1}{2}{\nu }^T{\nabla }^2{f}_i\left({\tilde{x}}_1\right)\nu \le \delta ,i=1,...,k,\hfill \\ \hfill & h_i\left({\tilde{x}}_1\right)+\nabla h_i{\left({\tilde{x}}_1\right)}^T\nu =0,i=1,...,p.\hfill \\ \hfill & g_i\left({\tilde{x}}_1\right)+\nabla g_i{\left({\tilde{x}}_1\right)}^T\nu =0,i\in I_c\left(ϵ\right).\hfill \end{array}$$

(24)

Note that problem $\left(24\right)$ is identical to problem $\left(23\right)$ except that $|I_c\left(ϵ\right)|$ inequalities are treated as equalities at ${\tilde{x}}_1$. In particular, we propose to add an index i to $I_c\left(ϵ\right)$ if

(a)

$g_i\left({\tilde{x}}_1\right)>ϵ$, i.e., if ${\tilde{x}}_1$ significantly violates the constraint $g_i$; or

(b)

$g_i\left({\tilde{x}}_1\right)\in (-ϵ,ϵ)$ and $\nabla g{\left({\tilde{x}}_1\right)}^T{\tilde{\nu }}_1>0$, i.e., if $x_i$ is either active but $g_i$ nearly active at $x_i$ or if $x_i$ already slightly violates $g_i$ and a step into direction ${\tilde{\nu }}_1$ would lead to (further) violation of this constraint, indicated by $\nabla g{\left({\tilde{x}}_1\right)}^T{\tilde{\nu }}_1>0$.

Algorithm 1 shows the pseudo code to build the index set $I_c\left(ϵ\right)$ at a predictor point ${\tilde{x}}_i$. Given the Newton direction, the Newton step can then be performed via using the Armijo rule described in [51], as done in our computations. The set $I_c\left(ϵ\right)$ is only computed once, it and remains fixed during the Newton iteration in the corrector step.

Algorithm 2 shows the pseudo code of one predictor–corrector step of the PT for general (equality and inequality constrained) MOPs. For bi-objective problems, $\mu$ can be chosen as in (16) leading either to a “left up” or “right down” movement, as discussed above. The algorithm has to be stopped if $\alpha$ is either close enough to ${(1,0)}^T$ or ${(0,1)}^T$, depending of course on the chosen search direction. For $k>2$, one can use the box partition in objective space as described in [2] in order to mark the regions of the Pareto front that have already been “covered” during the run of the algorithm.

For the realization of the predictor–corrector step several linear systems of equations have to be solved, the largest one being (20). The cost is hence $O\left({(n+p+s)}^3\right)$ in terms of flops and $O\left({(n+p+s)}^2\right)$ in terms of storage. Further, for the corrector step the SOP (7) has to be solved that contains $k+p+s$ decision variables. For the computation of the Newton direction, the SOPs (23) and (24) have to be solved for the first Newton iteration that contains both $n+1$ decision variables. For further Newton iterations, only SOP (24) has to be solved since the index set $I_c\left(ϵ\right)$ remains fixed within a corrector step. Finally, note that, if the method is realized as described above, the Hessians of all individual objectives have to be computed at each candidate solution (including at each Newton iteration). Using ideas from quasi-Newton methods, one can approximate the Hessians so that only gradient information is needed at each candidate solution, as described in [2].

Algorithm 1 Build $I_c\left(ϵ\right)$

Require:${\tilde{x}}_i$: predictor, $\widetilde{{\nu }_i}$: corrector direction for $\left(23\right)$, $ϵ>0$: tolerance Ensure:$I_c\left(ϵ\right)$: index set  1: $I:=\varnothing$  2: for$i=1,...,m$do  3:    if$g_i\left({\tilde{x}}_i\right)>ϵ$then  4:        $I:=I\cup i$  5:    else if$g_i\left({\tilde{x}}_i\right)\in (-ϵ,ϵ)\wedge \nabla g{\left({\tilde{x}}_i\right)}^T{\tilde{\nu }}_i>0$then  6:        $I:=I\cup i$  7:    end if  8: end for  9: Return$I_c\left(ϵ\right)$

Algorithm 2 Predictor–corrector step of the Pareto Tracer for general MOPs

Require:$x_i$: current candidate solution, $\tau >0$: desired distance in objective space, $ϵ>0$: tolerance Ensure:$x_{i+1}$: new candidate solution  1: Compute ${\alpha }^{\left(i\right)}\in {\mathbb{R}}^k$ via solving (7)  2: Compute $W_{\alpha }^{\left(i\right)}$ as in (19)  3: Compute A as in (10)  4: Select ${\mu }^{\left(i\right)}$ as in (16) or via (15) and (21)  5: Compute ${\nu }_{\mu }^{\left(i\right)}$ via solving (20)  6: $t^{\left(i\right)}:=\frac{\tau }{\parallel J{\nu }_{\mu }^{\left(i\right)}{\parallel }_2}$  7: ${\tilde{x}}_{i+1}:=x_i+t^{\left(i\right)}{\nu }_{\mu }^{\left(i\right)}$  8: Compute $x_{i+1}$ via a Newton method starting at ${\tilde{x}}_{i+1}$. For the Newton direction use the solution of (24).  9: Return $x_{i+1}$

As a demonstration example, we consider the problem

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{c}{f}_1\left(x\right)={(x_1+3)}^2+{(x_2-2)}^2,\hfill \\ f_2\left(x\right)=x_1^2+{(x_2+3)}^2,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& g_1\left(x\right)={(x_1+1)}^2+x_2^2\le 2^2,\hfill \\ \hfill & g_2\left(x\right)={(x_1+2)}^2+{(x_2+2)}^2\le 2^2.\hfill \end{array}$$

(25)

Figure 1a shows the Pareto set of the above problem where the two inequalities have been left out (i.e., the line segment connecting ${(-3,2)}^T$ and ${(0,-3)}^T$), the sets $g_i\left(x\right)=0$, $i=1,2$, as well as the Pareto set of this problem which is indeed the result of the PT. As starting point, we chose a point which significantly violates both constraints (and, hence, $|I_c\left(ϵ\right)|=2$ for $ϵ=1e-4$). An application of the above-described Newton method leads to the point on the Pareto set with the smallest $x_1$-value, which is in fact the initial point for the PT. During the run of PT, first only $g_2$ is “active” in the corrector step (i.e., $I_c\left(ϵ\right)=\left\{2\right\}$), later none of the constraints (in the intersection of the Pareto fronts of the constrained and the unconstrained MOP), and finally only $g_1$.

4. Numerical Results

In this section, we further demonstrate the behavior of the PT on five benchmark problems that contain inequality constraints. For all problems, we used the quasi-Newton variant of PT that only required function and Jacobian information (and no Hessians). To compare the results, we also show the respective results obtained by the normal boundary intersection (NBI, [16]), the $ϵ$-constraint method [5], and the multi-objective evolutionary algorithm NSGA-II. For NBI and the $ϵ$-constraint method, we used the code that is available at [52], and for NSGA-II the implementation of PlatEMO [53]. Regrettably, no comparison to a multi-objective continuation method can be presented since none of the respective codes are publicly available. For a comparison of the PT and the method of Hillermeier on box and equality constrained MOPs, we refer to [2]. We chose also to include a comparison to the famous NSGA-II since it is widely used and state-of-the-art for two- and three-objective problems as we consider here. We stress that the comparisons only show (on the first four test problems) that PT outperforms NSGA-II on these particular cases where the Pareto front consists of one connected component. For highly multi-modal functions where the Pareto set/front falls into several connected components, NSGA-II will certainly outperform the (standalone) PT. A fair comparison can only be obtained when integrating PT into a global heuristic (as, e.g., done in [41]). This is certainly an interesting task, however, beyond the scope of this work.

To compare the results, we compare the total number of function evaluations used for each algorithm on each problem. For this, each Jacobian call is counted as four function calls assuming that the derivative is obtained via automatic differentiation [54]. To measure the quality of the approximations, we used the averaged Hausdorff distance ${\mathsf{\Delta }}_2$ [55,56,57]. Since NSGA-II has stochastic components, we applied this algorithm for each problem 10 times and present the median result (measured by ${\mathsf{\Delta }}_2$).

4.1. Binh and Korn

Our first test example is a modification of the box-constrained BOP from Binh and Korn [58], where we add two inequality constraints as follows:

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{c}{f}_1\left(x\right)=4x_1^2+4x_2^2,\hfill \\ f_2\left(x\right)={(x_1-5)}^2+{(x_2-5)}^2,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {(x_1-2)}^2+{(x_2-1)}^2\le 2.3^2,\hfill \\ \hfill & {(x_1-3)}^2+{(x_2-3)}^2\ge 1.5^2,\hfill \\ \hfill & 0\le x_1\le 5,\hfill \\ \hfill & 0\le x_2\le 3.\hfill \end{array}$$

(26)

Table 1. Parameters used by NSGA-II for MOP (26).

Population Size	100
Number of generations	20
Probability of crossover	0.9
Probability of mutation	0.5

shows the design parameters that have been used by NSGA-II for this problem,

Table 2. Computational efforts and approximation quality of the algorithms for MOP (26).

	PT	NBI	$\mathit{ϵ}$-Constr.	NSGA-II
Solutions	52	52	52	100
Function Evaluations	151	427	336	2000
Jacobian Evaluations	133	425	336	-
Hessian Evaluations	-	373	284	-
Total of Evaluations	683	8095	6224	2000
${\mathsf{\Delta }}_2$	0.6050	0.6025	0.9272	0.5626

shows the computational efforts and the obtained approximation quality for each algorithm, and Figure 2 and Figure 3 show the obtained Pareto set and front approximations, respectively. For PT, we chose $\tau =0.6$ leading to 52 solutions along the Pareto set/front in 4.48 s (the computations have been done on a Ubuntu 20.04.1 LTS system with an Intel Core i7-855OU 1.80 GHz x 8 CPU and 12 GB of RAM). We then applied NBI and the $ϵ$-constraint model using this number of sub-problems. For NSGA-II, we took the population size 100, which is a standard value for this algorithm. The results show nearly perfect Pareto front approximations (at least from the practical point of view) for all algorithms, which is also reflected by the low ${\mathsf{\Delta }}_2$ values that are very close to the optimal value of $0.6$ (at least for PT, defined by $\tau$). In terms of function evaluations, PT clearly wins over NBI and the $ϵ$-constraint method. A comparison to NSGA-II is not possible due to the choice of the population size.

4.2. Chakong and Haimes

Next, we consider the bi-objective problem of Chankong and Haimes [59], which contains next to the box constraints one linear and one nonlinear inequality.

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{c}{f}_1\left(x\right)=2+{(x_1-2)}^2+{(x_2-1)}^2,\hfill \\ f_2\left(x\right)=9x_1-{(x_2-1)}^2,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& x_1^2+x_2^2\le 225,\hfill \\ \hfill & x_1-3x_2+10\le 0,\hfill \\ \hfill \mathrm{with}& -20\le x_1,x_2\le 20.\hfill \end{array}$$

(27)

Table 3. Parameters used by NSGA-II for problem (27).

Population Size	100
Number of generations	30
Probability of crossover	0.9
probability of mutation	0.5

shows the parameter values used for the application of NSGA-II,

Table 4. Computational efforts and approximation qualities for problem (27).

	PT	NBI	$\mathit{ϵ}$-Constr.	NSGA-II
Solutions	80	80	80	100
Function Evaluations	540	678	578	3000
Jacobian Evaluations	499	678	578	-
Hessian Evaluations	-	598	498	-
Total of Evaluations	2536	12,958	10,858	3000
${\mathsf{\Delta }}_2$	1.1459	1.1457	1.2141	1.1871

the computational efforts and the approximation qualities, and Figure 4 and Figure 5 the obtained approximations. We used $\tau =1$ for PT, and proceeded as for the previous example for the other methods. The results are also similar to the previous example: all methods are capable of detecting a nearly perfect Pareto front approximation, and the overall cost is significantly less for PT, in 5.96 s.

4.3. Tamaki

Next, we considered a MOP with three objectives (28):

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{cc}{f}_1\left(x\right)\hfill & =x_1,\hfill \\ f_2\left(x\right)\hfill & =x_2,\hfill \\ f_3\left(x\right)\hfill & =x_3,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& x_1^2+x_2^2+x_3^2\ge 1,\hfill \\ \hfill & 0\le x_1,x_2,x_3\le 4.\hfill \end{array}$$

(28)

Both the Pareto set and front for this problem are a part of the unit sphere.

Table 5. Parameters used by NSGA-II for MOP (28).

Population Size	300
Number of generations	150
Probability of crossover	0.9
Probability of mutations	0.5

shows the design parameters for NSGA-II,

Table 6. Computational efforts and approximation qualities for problem (28).

	PT	NBI	$\mathit{ϵ}$-Constr.	NSGA-II
Solutions	305	112	N/A	300
Function Evaluations	2498	3758	N/A	450,000
Jacobian Evaluations	1101	3758	N/A	-
Hessian Evaluations	-	3293	N/A	-
Total of Evaluations	6902	91,236	N/A	450,000
${\mathsf{\Delta }}_2$	0.0380	0.6353	N/A	0.0390

shows the computational effort and the approximation quality for each algorithm, and Figure 6 shows the Pareto front approximations (the respective Pareto set approximations will look identically, albeit in x-space). For this problem, $\tau =0.05$ was used. The implementation of the $ϵ$-constrained method did not yield a result. On the Tamaki problem, PT performs better than the other algorithms both in approximation quality and in the overall computational cost.

4.4. BCS

We next considered a second three-objective problem that contains next to one inequality also a linear equality constraint:

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{c}{f}_1\left(x\right)={(x_1+3)}^2+{(x_2+3)}^2+{(x_3+3)}^2,\hfill \\ f_2\left(x\right)={(x_1-9)}^2+{(x_2+5)}^2+{(x_3+5)}^2,\hfill \\ f_3\left(x\right)={(x_1-5)}^2+{(x_2-8)}^2+x_3^2,\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& x_1-2x_2-3x_3=0,\hfill \\ \hfill & sin\left(2x_1\right)-x_2\le 0.\hfill \end{array}$$

(29)

Table 7. Parameters used by NSGA-II for MOP (29).

Population Size	100
Number of generations	500
Probability of crossover	0.9
Probability of mutations	0.5

presents the design parameters used by NSGA-II,

Table 8. Computation efforts for the proposed test problem (29).

	PT	NBI	$\mathit{ϵ}$-Constr.	NSGA-II
Solutions	378	0	N/A	4
Function Evaluations	1923	2290	N/A	50,000
Jacobian Evaluations	756	1641	N/A	-
Hessian Evaluations	-	1431	N/A	-
Total of Evaluations	4947	40,336	N/A	50,000
${\mathsf{\Delta }}_2$	2.0658	-	N/A	61.7685

shows the computational effort and the approximation quality for each algorithm, and Figure 7 and Figure 8 present the Pareto front approximation of PT (using $\tau =2$), which took 16.79 s. For this example, none of the other methods were able to yield feasible solutions, where we counted a solution x to be feasible if $|x_1-2x_2-3x_3|<1e-4$ and $sin\left(2x_1\right)-x_2\le 1e-4$.

4.5. Osykzka and Kundu

As last example, we considered the bi-objective problem of Osykzka and Kundu [60], which has six decision variables and contains six inequality constraints in addition to the box constraints:

$$\begin{array}{cc}\hfill min& \left\{\begin{array}{cc}{f}_1\left(x\right)\hfill & =-25{(x_1-2)}^2-{(x_2-2)}^2-{(x_3-1)}^2-{(x_4-4)}^2-{(x_5-1)}^2\hfill \\ f_2\left(x\right)\hfill & {\displaystyle =\sum _{i=1}^6{x}_i^2}\hfill \end{array}\right.\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& x_1+x_2-2\ge 0\hfill \\ \hfill & 6-x_1-x_2\ge 0\hfill \\ \hfill & 2-x_2+x_1\ge 0\hfill \\ \hfill & 2-x_1+3x_2\ge 0\hfill \\ \hfill & 4-{(x_3-3)}^2-x_4\ge 0\hfill \\ \hfill & {(x_5-3)}^2+x_6-4\ge 0\hfill \\ \hfill & 0\le x_1,x_2,x_6\le 10\hfill \\ \hfill & 1\le x_3,x_5\le 5\hfill \\ \hfill & 0\le x_4\le 6\hfill \end{array}$$

(30)

While the Pareto front of this problem is connected, its Pareto set consists of three different connected components. Hence, PT is not able to compute an approximation of the entire Pareto front with only one starting point. Figure 9a shows the result of PE for $\tau =2$ using the three starting points

$$\begin{array}{cc}\hfill x_{0,1}& ={(0.60,1.50,1.0,0.00,1.00,0.04)}^T,\hfill \\ \hfill x_{0,2}& ={(0.00,2.00,2.20,0.00,1.00,0.00)}^T,\hfill \\ \hfill x_{0,3}& ={(5.00,1.00,5.00,0.00,1.00,0.01)}^T.\hfill \end{array}$$

(31)

The computational time to obtain this result was 12.98 s. Figure 9b shows a numerical result of NSGA-II using the design parameters shown in

Table 9. Parameters used by NSGA-II for MOP (30).

Population Size	435
Number of generations	50
Probability of crossover	0.9
Probability of mutations	0.5

. The obtained solutions “under” the Pareto front can be explained by the tolerance of $1\times {10}^{-4}$ that was used to measure feasibility (while $1\times {10}^{-8}$ was used for PT).

Table 10. Computational efforts and approximation qualities for problem (30).

	PT	NSGA-II
Solutions	435	428
Function Evaluations	2051	20,000
Jacobian Evaluations	850	-
Hessian Evaluations	-	-
Total of Evaluations	5451	20,000
${\mathsf{\Delta }}_2$	1.801	2.4819

shows the computational effort for both methods. Needless to say, this represents by no means a comparison of the two methods. Instead, this should be rather seen as a motivation to hybridize PT with a global search strategy in order to obtain a fast and reliable multi-objective solver, which we leave for future studies.

5. Conclusions and Future Work

In this paper, we extend the multi-objective continuation method Pareto Tracer (PT) for the treatment of general inequality constraints. To this end, the predictor–corrector step is modified as follows: in the predictor, all nearly active inequalities are treated as equalities. In the following corrector step, the main challenge is to identify the inequalities for which the predictor solution is either nearly active or slightly violates the constraint that has to be considered, namely the equality constraint in the Newton method, and this is done in a bootstrap manner. We formulate the resulting algorithm and show some numerical results on several benchmark problems, indicating that it can reliably handle inequality (and equality) constrained MOPs. We further present comparisons to some other numerical methods. The results show that the extended PT can indeed reliably handle general MOPs (and in particular general inequalities). However, the method is—by construction—of local nature and restricted to the connected component of the solution set for which one initial solution is available. One interesting task is certainly to hybridize PT with a global solver such as a multi-objective evolutionary algorithm and to compare the resulting hybrid against other methods with respect to their ability to compute the entire global Pareto set/front of a given MOP. This is beyond the scope of this work and has been left for future work.