Algorithms for Generating Pareto Fronts of Multi-objective Integer and Mixed-Integer Programming Problems

TL;DR

This paper adapts algorithms originally for continuous problems to approximate Pareto fronts in multi-objective integer and mixed-integer programming, demonstrating effectiveness on complex problems including rocket injector design.

Contribution

It introduces adapted algorithms for multi-objective integer and mixed-integer problems and analyzes scalarization techniques to improve Pareto front approximation.

Findings

Algorithms successfully approximate Pareto fronts for integer problems.

Effective handling of multi-objective mixed-integer problems demonstrated.

New four-objective algorithm applied to rocket injector design.

Abstract

Multi-objective integer or mixed-integer programming problems typically have disconnected feasible domains, making the task of constructing an approximation of the Pareto front challenging. The present paper shows that certain algorithms which were originally devised for continuous problems can be successfully adapted to approximate the Pareto front for integer, and mixed-integer, multi-objective optimization problems. Relationships amongst various scalarization techniques are established to motivate the choice of a particular scalarization in these algorithms. The proposed algorithms are tested by means of two-, three- and four-objective integer and mixed-integer problems, and comparisons are made. In particular, a new four-objective algorithm is used to solve a rocket injector design problem with a discrete variable, which is a challenging mixed-integer programming problem.

Tables

Table 1. Table 1: The algorithms implemented in the paper.

	Number of	Grid
Algorithm	objectives	generation	Scalarization	Reference
1	2	CHIM	Weighted-constraint	[7]
2	2	CHIM	Pascoletti–Serafini	[7]
3	3	CHIM	Weighted-constraint	[8]
4	3	CHIM	Pascoletti–Serafini	[8]
5	3	SBG	Weighted-constraint	[8]
6	3	SBG	Pascoletti–Serafini	[8, 34]
7	4	SBG	Weighted-constraint	New

Table 2. Table 2: Test problem 1 – Numerical performance of Algorithms 1 and 2.

		Number of	Number of
	CPU time	subproblems	weak Pareto
Algorithm	[sec]	attempted	points generated
1	14	$11\times 2=22$	9 (all)
2	13	14	6

Table 3. Table 3: Test problem 2 – Numerical performance of Algorithms 3 and 4.

		Number of	Number of
	CPU time	subproblems	weak Pareto
Algorithm	[sec]	attempted	points generated
3	147	$44\times 3=132$	19 (all)
4	150	135	7

Table 4. Table 4: Test problem 3 – Numerical performance of Algorithms 3–6.

		Number of	Number of
	CPU time	subproblems	weak Pareto
Algorithm	[sec]	attempted	points generated
3	202	$230\times 3=690$	58
4	211	495	33
5	$38+8+82=128$	$45+91+273=409$	60 (all)
6	$32+14+84=130$	$60+170+170=400$	33
6	$184+1400+5530=7114$	$360+3081+3081=6522$	44

Equations

$$\left\{\begin{array}[]{rl}\displaystyle\min&\ \ f(x):=\left[f_{1}(x),\ldots,f_{\ell}(x)\right]\\[5.69054pt] \mbox{s.t.}&\ \ x\in X:=\left\{x\in\mathbb{Z}^{n_{1}}\times\mathbb{R}^{n_{2}}\mid g_{j}(x)\leq 0,\;\;j=1,\ldots,m\right\},\end{array}\right.$$

$$\underset{i=1,\dots ,\ell }{\min }\left\{\underset{x\in X}{\min }{f}_i(x)\right\}>0.$$

$$W:=\left\{w\in {\mathbb{R}}^{\ell }\mid w_i>0,\sum _{i=1}^{\ell }{w}_i=1\right\}.$$

$$\underset{x\in X}{\min }{f}_i(x).$$

$$\underset{x\in X}{\min }\sum _{i=1}^{\ell }{w}_i{f}_i(x).$$

$$Sol(P_w):=\{x\in X\mid x\text{mbox}solves(P_w)\}.$$

$$\left\{\begin{array}[]{cl}\displaystyle\min_{x\in X}&\ \ f_{k}(x)\\[5.69054pt] \mbox{s.t.}&\ \ f_{i}(x)\leq\epsilon_{i}\,,\ \ \forall\;i=1,\ldots,\ell\,,\;\;i\neq k,\end{array}\right.$$

$$Sol(P^k(ϵ)):=\{x\in X\mid x\text{mbox}solves(P^k(ϵ))\}.$$

$$\left\{\begin{array}[]{cl}\displaystyle\min_{(\alpha,x)\in\mathbb{R}\times X}&\ \ \alpha\\[5.69054pt] \mbox{s.t.}&\ \ w_{i}\,(f_{i}(x)-u_{i})\leq\alpha\,,\ \ \forall\;i=1,\ldots,\ell\,,\end{array}\right.$$

$$\left\{\begin{array}[]{cl}\displaystyle\min_{x\in X}&\ \ f_{k}(x)\\[5.69054pt] \mbox{s.t.}&\ \ w_{i}\,f_{i}(x)\leq w_{k}\,f_{k}(x)\,,\quad\forall\;i=1,\ldots,\ell,\;\;i\neq k\,,\end{array}\right.$$

$$Sol(P_w^k):=\left\{x\in X\mid x\text{mbox}solves(P_w^k)\right\}.$$

$$\bigcap _{k=1}^{l}Sol(P_{w^{\prime }}^k)\subset WE(P)\subset \bigcup _{w\in W}\left[\bigcap _{k=1}^{l}Sol(P_w^k)\right].$$

$$w_i:=\frac{1/f_i(x)}{\sum _{j=1}^{\ell }1/f_j(x)},i=1,\dots ,\ell .$$

$$x\in \bigcap _{k=1}^{\ell }Sol(P_w^k),$$

$$f_r({\bar{x}}_r)\ge f_r({\bar{x}}_k).$$

$$w_i:=\frac{1/f_i(\hat{x})}{\sum _{j=1}^{\ell }1/f_j(\hat{x})},i=1,\dots ,\ell .$$

$$w_k{f}_k(\hat{x})=w_i{f}_i(\hat{x}),\forall i=1,2,...,\ell ,i\ne k.$$

$$\ \left\{\begin{array}[]{cl}\displaystyle\min_{x\in X}&\ \ w_{k}f_{k}(x),\\[11.38109pt] \mbox{s.t.}&\ \ w_{i}f_{i}(x)\leq w_{k}f_{k}(\hat{x}),\ \ i=1,2,...,\ell\,,\ i\neq k.\end{array}\right.$$

$$Sol(P_{\hat{x}}^k):=\{x\in X\mid x\text{mbox}solves(P_{\hat{x}}^k)\}.$$

$$w_i{f}_i(x)\le w_i{\widehat{ϵ}}_i=w_i{f}_i(\hat{x})=w_k{f}_k(\hat{x}),$$

$$\sum _{i=1}^{\ell }{w}_i[f_i(x)-f_i(\bar{x})]\ge 0,\forall x\in X.$$

$${\alpha }_j{f}_j(x)<{\alpha }_j{f}_j(\bar{x}),$$

$${\alpha }_i{f}_i(x)\le {\alpha }_j{f}_j(\bar{x})={\alpha }_i{f}_i(\bar{x}),i=1,2,...,\ell \text{mbox}andi\ne j,$$

$$f_j(x)-f_j(\bar{x})<0,$$

$$f_i(x)-f_i(\bar{x})\le 0,i\ne j,$$

$$\sum _{i=1}^{\ell }{w}_i(f_i(x)-f_i(\bar{x}))<0,$$

$$f_k(\hat{x})<f_k(\bar{x})\text{mbox}forsomek.$$

$$f_k(\bar{x})\le f_k(x),\forall x\in X\text{mbox}suchthat{f}_i(x)\le {\widetilde{ϵ}}_i,i\ne k.$$

$$f_k(\bar{x})\le f_k(\hat{x}),$$

$$\begin{array}[]{cl}\rm\min&\ (x_{1},\ x_{2})\\[5.69054pt] \mbox{s.t.}&(x_{1}-4)^{2}+(x_{2}-4)^{2}-16\leq 0\,,\\[2.84526pt] &0\leq x_{1},x_{2}\leq 4\,,\;\;\mbox{and}\;\;x_{1},x_{2}\mbox{ are integers}.\end{array}$$

Full text

Key words: Multi-objective optimization, Integer programming, Mixed-integer programming, Scalarization, Pareto front, Efficient set, Numerical methods.

AMS subject classifications. 90C26, 90C29, 90C30, 90C56.

1 Introduction

Multi-objective optimization is concerned with simultaneous minimization of several conflicting objective functions. The aim is to obtain/approximate, and view, the set of trade-off, or compromise, solutions. This set of trade-off solutions is referred to, in the current article, as the (weak) Pareto front, or the (weak) efficient set, of the problem—a more precise definition is to be provided in Section 2. Once the set of Pareto points is at hand, the decision maker conveniently selects a suitable solution from this set, usually based on an additional criterion. These kinds of problems arise in many multidisciplinary applications, for example, in health [52, 9, 61], engineering [30, 50], mining [59, 26], and finance [45, 49].

If all variables of a multi-objective optimization problem belong to a continuous set, then the problem is called a multi-objective continuous optimization problem. If, on the other hand, the variables belong to a subset of integers, then the problem is a multi-objective integer programming one. When the model involves both integer and continuous variables the problem is referred to as a multi-objective mixed-integer programming problem. The latter case arises in many real-life applications, for example, the knapsack [28, 10] and shortest-path problems [13, 54]. Other types of mixed-integer multi-objective problems can be found in [43].

It is desirable that algorithms for solving multi-objective integer and mixed-integer programming problems will have the following major practical attributes: (i) construct a faithful approximation of the Pareto front and (ii) generate this approximation in a reasonable amount of computational time. Attribute (ii) is particularly valuable, since, for example, solving each scalarized problem for a mixed-integer programming problem can be very costly.

Indeed it is in general not realistic to compute too many points in the Pareto front, especially if the Pareto front contains connected segments (e.g. segments of curves in the case of two objectives and segments of surfaces in the case of three objectives). Therefore we aim to construct an approximation of the Pareto front in the following sense. We consider some discretization grid, or a finite partition, of the space of all weights for a parametric scalarization of the multi-objective problem.

In Sections 3.1–3.3, relationships amongst various scalarization techniques are presented. Namely, it is shown that every point in the Pareto front that can be generated by the weighted-sum [19, 60], $k$th-objective $\epsilon$-constraint [11] or feasible-value-constraint scalarization [8] technique, can also be generated by the weighted-constraint scalarization technique introduced in [8] and implemented in various algorithms in [8, 7]. On the other hand, a solution of the weighted-constraint scalarization may not be a solution of the other three scalarization techniques. Therefore, in the numerical experiments, the weighted-sum, $k$th objective $\epsilon$-constraint and feasible-value-constraint scalarization techniques are not incorporated into the algorithms that are used to approximate the Pareto front.

As a result, the scalarization techniques we use assure that the solution of the scalarized problem with each weight in the grid yields a Pareto solution, and the collection of all these solutions constitutes a finite and faithful subset of the Pareto front. This finite subset is an approximation generated by the weak efficient points; in other words, all the points in this finite set are non-dominated by each other.

The phrase, approximation of the Pareto front, is used in the sense that, as the discretization is made finer and finer, a larger subset of the true Pareto front is generated. The same notion of approximation was used in the authors’ previous studies in [7, 8].

There is an abundance of practical algorithms for constructing an approximation of the Pareto front of continuous-variable problems, especially those with two or three objective functions [7, 8, 34, 15, 17, 18, 27, 38, 44, 57, 46]. A successful algorithm in the case of four objectives has recently been proposed and numerically studied by the authors of the current paper in [8].

As for integer and mixed-integer programming problems, the question of generating a faithful approximation of the Pareto front is still a developing area for problems with three or more objectives. Multi-objective integer problems (MOIP), where all variables are integer, have been studied and effective algorithms proposed in earlier and recent literature, see e.g. [32, 14, 4, 16]. The present paper focuses more on the mixed-integer case, acknowledging the fact that mixed-integer programming problems are in general far more challenging than the integer ones. For problems which are posed in convex form, the reader is referred to the branch-and-bound approach presented in the recent article [16]. Those approaches that do work with general three objectives, do not attempt to construct the whole front—see [2, 3, 31, 36, 37, 40, 41, 42, 35]. Hence there is a need for new techniques that approximate Pareto fronts for the mixed-integer case. Given the fact that for these type of problems there is no certificate of optimality, to solve them numerically, one needs a suitable solver for the subproblems that guarantees to yield a globally optimal solution. To the best of the authors’ knowledge, such a solver is not (yet) available for mixed-integer programming problems.

Presumably, one of the main reasons for which there are almost no studied examples of three- and four-objective mixed-integer optimization problems in the literature is that the feasible set in the presence of integer variables is often disconnected. In this situation, specialized algorithms are needed to construct the Pareto front. Specialized scalarization techniques for problems with continuous variables and disconnected domains were first studied by the authors in [44, 7] and then in [8]—also see [38, Sections 2.4 and 3.2.4].

Approximation of the Pareto front indeed relies on the success of the scalarization technique employed to transform the problem into a single variable one. If the multi-objective problem is continuous, and the domain is connected, then the Pascoletti–Serafini scalarization [39, 18], or the constrained-Tchebycheff scalarization [17], can be successfully used to construct the Pareto front. If the domain is disconnected, however, the Pascoletti–Serafini scalarization, as well as most of the other available scalarization techniques, are likely to fail, as illustrated in [8, 7]. In such cases, implementation of certain other scalarization techniques, in particular the weighted-constraint scalarization technique and the associated algorithms introduced in [8, 7], have been shown to perform well in constructing the Pareto front.

The algorithms in [8, 7] have been demonstrated to be particularly successful when applied to three- and four-objective continuous optimization problems. However, the same algorithms have never been tested before on multi-objective integer or mixed-integer optimization problems. Therefore it is natural to extend and implement the algorithms in [8, 7] for mixed-integer problems. If a problem is linear then the solver PolySCIP developed by Schenker [48] or the search based algorithms introduced in [5], can possibly be used. Hence, there is a need to devise algorithms for nonlinear and nonconvex mixed-integer optimization problems, in particular for problems with more than three objective functions. Note that an integer or mixed-integer problem is called nonconvex when the objective function or the form of the constraints is not convex.

The current paper shows that Algorithms 1 [7], 3 and 5 [8] are particularly successful in solving multi-objective integer or mixed-integer programming problems, compared with Algorithms 2 [7], 4 [8] and 6 [34, 8]. The latter group of algorithms use the conventional Pascoletti–Serafini scalarization, which is effective in the case of simply connected Pareto fronts. Moreover, the new Algorithm 7 is an extension of Algorithm 5 to the case of four objective functions. Algorithm 7 itself and its application to a rocket injector design problem [8, 21, 55] in Section 5.5 is the first of its kind, as, to the knowledge of the authors, there does not exist a multi-objective mixed-integer optimization algorithm applicable to problems with three or more nonlinear objective functions. See Table 1 and Section 4.2 for further details.

The main contributions of the present paper are summarized as follows.

• •

A theoretical comparison study of the four scalarization techniques mentioned above is carried out so as to motivate the choice of the particular techniques employed in the numerical implementation of the algorithms.

• •

Algorithm 5 given in [8] is extended to the case of four-objective mixed-integer programming problems, and implemented for a challenging test problem.

• •

Algorithm 1 given in [7]—also see [38, Section 3.2.4] and Algorithm 8.1 given in [34] are tested and compared for three-objective integer programming problems. Note that these algorithms were previously tested in [7] and [34] only on continuous optimization problems.

The layout of the paper is as follows. The problem description along with some preliminaries are presented in Section 2. Five scalarization approaches, namely the Pascoletti–Serafini, weighted-sum, $k$th-objective $\epsilon$-constraint, feasible-value-constraint, and weighted-constraint scalarization, are recalled in Section 3. In this section, relationships among several of these scalarization techniques are established. In Section 4, algorithms for two and three objectives are described, and a new algorithm for four objectives is introduced. Numerical experiments with these algorithms are given in Section 5. In particular, Section 5.5 presents the results for a new and challenging reformulation of the rocket injector design problem involving four objective functions. Conclusion and discussion is given in Section 6. Detailed descriptions of Algorithms 3 and 7 are included in the Appendix.

2 Preliminaries

In this section, we collect the relevant notions, definitions, and concepts that are used in our study. These are standard notations and tools for multi-objective optimization, and we follow here the classical notation in the literature—see for example, [11, 33, 58]. Let $\mathbb{R}^{n}$ be the $n$-dimensional Euclidean space. Define $\mathbb{R_{+}}$ to be the set of all nonnegative real numbers and $\mathbb{R_{++}}$ the set of positive numbers. Given the positive integers $\ell,m,n_{1}$ and $n_{2}$, we consider the following multi-objective optimization Problem $(P)$.

[TABLE]

where $\mathbb{Z}$ is set of all integers, and $f_{i}:\mathbb{Z}^{n_{1}}\times\mathbb{R}^{n_{2}}\to\mathbb{R}$, $i=1,\ldots,\ell$, and $g_{j}:\mathbb{Z}^{n_{1}}\times\mathbb{R}^{n_{2}}\to\mathbb{R}$, $j=1,\ldots,m$. We assume that the functions $f_{i}$ are bounded below on the constraint set $X$, and the lower bounds of the functions $f_{i}$ are known. Therefore, we can impose that

[TABLE]

We will use the following set of positive weights in our analysis.

[TABLE]

We recall next two types of solutions of $(P)$. The more restrictive type of solution is the so-called efficient point [58] or Pareto point [33], and a less restrictive concept is the one of a weak efficient point. We also provide below the concepts of ideal and utopia vectors and individual minima.

Definition 2.1

(See [11] and [33])

• (a)

A point $\bar{x}\in X$ is said to be efficient for Problem $(P)$ iff there is no $x\in X$, such that $f_{i}(x)\leq f_{i}(\bar{x})$, $\forall\,i\in\{1,\ldots,\ell\}$, and $f_{j}(x)<f_{j}(\bar{x})$, for some $j\in\{1,\ldots,\ell\}$. Let, $E(P)$ be the set of efficient points of Problem $(P)$.

• (b)

A point $\bar{x}\in X$ is said to be weak efficient for Problem $(P)$ iff there is no $x\in X$ such that $f_{i}(x)<f_{i}(\bar{x})$, $\forall\,i\in\{1,\ldots,\ell\}$. Let, $WE(P)$ be the set of weak efficient points of Problem $(P)$. We define the Pareto front of Problem $(P)$ as the image of $WE(P)$ under $f$.

• (c)

Suppose that $\bar{x}_{f_{i}}$ is a minimizer of $f_{i}$, $i=1,\ldots,\ell$, over the set $X$. That is to say, $\bar{x}_{f_{i}}$ solves the optimization problem

[TABLE]

Then the vector $\bar{x}_{f}=[\bar{x}_{f_{1}},\ldots,\bar{x}_{f_{\ell}}]$ is called an ideal vector. We denote the vector of individual minima (IM) by $f(\bar{x}_{f_{i}})=[f_{1}(\bar{x}_{f_{1}}),f_{2}(\bar{x}_{f_{2}}),\ldots,f_{\ell}(\bar{x}_{f_{\ell}})]$.

• (d)

A utopia vector $u:=(u_{1},\ldots,u_{\ell})\in\mathbb{R}^{\ell}$ associated with Problem $(P)$ is defined by $u_{i}:=f_{i}(\bar{x}_{f_{i}})-\varepsilon_{i}$, where, for all $i=1,\ldots,\ell$, $\varepsilon_{i}>0$.

For the problem to be truly multi-objective, ideal vectors cannot be feasible. It is worth noting that $E(P)$ $\subseteq$ $WE(P)$, but the opposite inclusions in general do not hold.

At first glance, weak efficient solutions may not appear to be convenient from a practical point of view, since they could (in theory) be improved. On the other hand, some Pareto fronts are extremely difficult to approximate (as in the case of the rocket injector design problem that we study in Section 5.5), and practitioners may have to look for weak efficient solutions which are not efficient as well as the efficient ones. In other words, in situations when finding an efficient solution is not an easy task, weak efficient solutions need to be welcome.

Weak efficient points may also help eliminate (or filter out) some of the dominated points produced by the numerical procedure. In particular, they allow one to obtain necessary conditions for non-dominated solutions—see [23, 29].

Besides from these valuable theoretical aspects of weak efficient points, we refer to the works [58, 25, 51, 11], where the approximation of the whole Pareto Front (efficient and weak efficient points) has important applications.

Studies of existence of weak efficient points, such as the ones in [22, 24] demonstrate their relevance in multi-objective optimization. A concrete example of the usefulness of approximating the entire Pareto front (efficient and weak efficient points), is illustrated later on via one of our test problems in Remark 5.1.

A comment is in order regarding the different parts of the Pareto front. In what follows, we will be using the expressions interior and/or boundary of the Pareto front. We clarify next what we mean exactly by these latter two expressions, since they are not to be confused with the classical concepts of interior and boundary used in topology theory.

Since the Pareto front has no interior in $\mathbb{R}^{\ell}$ (otherwise we could find points in the front which are strictly dominated), it is homeomorphic to a subset of $\mathbb{R}^{\ell-1}$. Hence, when we refer to the boundary of the Pareto front, what we actually mean is the subset of the front (in the $\mathbb{R}^{\ell}$ space), which is homeomorphic to the topological boundary of its homeomorphic image in $\mathbb{R}^{\ell-1}$. Similarly, by interior of the Pareto front, we mean the subset of the front (in the $\mathbb{R}^{\ell}$ space), which is homeomorphic to the topological interior of its homeomorphic image in $\mathbb{R}^{\ell-1}$.

For instance, if $\ell=2$, the boundary of the Pareto front (if connected) consists of the points $(\min_{x\in X}f_{1}(x),\max_{x\in X}f_{2}(x))$ and $(\min_{x\in X}f_{2}(x),\max_{x\in X}f_{1}(x))$. If $\ell=3$ and the Pareto front in $\mathbb{R}^{3}$ is connected, then it is homeomorphic to the two-dimensional simplex. In this case, the boundary of the Pareto front is the union of regions homeomorphic to lower dimensional faces of the simplex. To avoid the use of involved terminology, we refer in what follows, to the boundary and interior of the Pareto front via this homeomorphism.

3 Scalarization Techniques and Their Relationships

In this section, we recall five scalarization techniques; namely, the

• (i)

weighted-sum scalarization [19],

• (ii)

$k$th-objective $\epsilon$-constraint scalarization [11],

• (iii)

Pascoletti–Serafini scalarization [39],

• (iv)

($k$th-objective) weighted-constraint scalarization [7],

• (v)

feasible-value-constraint scalarization [8].

The first three of these have been popular/classical approaches for problems with only continuous variables. The latter two were proposed by the authors of the current paper and tested only on problems with continuous variables. However, the weighted-constraint scalarization has been shown to be efficient for approximating Pareto points when the domain and/or front might be disconnected [7, 38]. Even though these techniques exist in the literature, for the convenience of the reader, we give here short descriptions. For more details, the reader is referred to [7, 8], as well as [38, Sections 2.4 and 3.2.4], and the references therein.

Weighted sum scalarization.  This type of scalarization was introduced by Gass and Saaty in [19] and Zadeh in [60], it is computationally cheap and easy to implement. However, it cannot generate efficient points lying in nonconvex111The word nonconvex is loosely used here to describe those sections of the front which are “dented”, or, “going inwards”. sections of the front. The idea behind this method is to minimize a weighted sum of all objectives. Namely, with fixed $w\in W$, the weighted-sum scalarization of Problem $(P)$ is given by

[TABLE]

Every solution of (Pw) is weak efficient (see [33]), and this fact is used in constructing an approximation of the Pareto front. For a fixed $w\in W$, define the solution set of Pw as

[TABLE]

The $k$th-objective $\epsilon$-constraint scalarization.  The method was introduced by Changkong and Haimes, and a comprehensive analysis can be found in [11]. Let $\epsilon:=(\epsilon_{1},\ldots,\epsilon_{\ell})\in\mathbb{R}^{\ell}$ be fixed. In this approach, each of the objectives (say, the function $f_{k}$), is minimized and the remaining $(\ell-1)$ objectives are constrained by upper bounds $\epsilon_{i}\in\mathbb{R}$, $i=1,...,\ell$, $i\neq k$. The $k$th-objective $\epsilon$-constraint scalarization for Problem $(P)$ is defined as

[TABLE]

where $k\in\{1,2,...,\ell\}.$ Unlike the weighted-sum approach, this method can generate the Pareto points located in nonconvex sections of the front. However, the upper bounds $\epsilon_{i}$ in the constraints have to be chosen carefully, otherwise the new feasible region might be empty. In order to avoid this situation, a suitable range of values for $\epsilon_{i}$ has to be known beforehand. In [33, p. 85–86], it is shown that if $\bar{x}$ solves ($P^{k}(\epsilon)$) for some $k$ then $\bar{x}\in WE(P)$. For $\epsilon\in\mathbb{R}^{\ell}$, we define

[TABLE]

Pascoletti–Serafini scalarization.  Introduced by Pascoletti and Serafini in [39], this approach is also referred to as the goal-attainment method (see [12, 17, 33, 39]). Given a fixed $w\in W$, the Pascoletti–Serafini scalarization is posed as the following problem.

[TABLE]

where $\alpha\in\mathbb{R}$ is a new variable and $u=(u_{1},\ldots,u_{\ell})$ is a utopia vector. It should be noted that Problem $(PS)$ can be reformulated as the (weighted) Tchebychev(-norm) scalarization [17]. Every solution of Problem $(PS)$ is a weak efficient point [18]. Hence, this technique is used for approximating connected Pareto fronts in [34]. A modified version of the Pascoletti–Serafini scalarization has been used in [15], again for generating connected Pareto fronts.

Weighted-constraint scalarization.  Introduced in [7] by the authors of the current paper, this scalarization has been numerically illustrated to be useful in finding weak efficient points in a disconnected Pareto front, in the presence of a disconnected feasible set or a disconnected domain. As in the $k$th objective $\epsilon$-constraint approach, each objective $f_{k}$, $k=1,\ldots,\ell$, is minimized separately; however, all of the objectives appear in the constraints of the scalarized problem. Namely, the following $\ell$ problems are solved.

[TABLE]

where $k=1,\ldots,\ell$. Each Problem $(P_{w}^{k})$, for $k=1,\ldots,\ell$, is referred to as a subproblem of the scalarization. For fixed $w\in W$ and $k\in\{1,\ldots,\ell\}$, define the solution set of a subproblem as

[TABLE]

Given $w^{\prime}\in W$, it is proved in [7] that

[TABLE]

When the leftmost expression is nonempty for some $w^{\prime}\in W$, the leftmost inclusion above can be used to obtain points in the Pareto front. If such a weight $w^{\prime}$ is not directly available then [7, Theorem 3.1 and Proposition 3.3] can be used to generate new points in the front. To make these statements precise, we recall next the relevant results. The first of these results ensures that the map from the space of weights $W$ to the Pareto set $WE(P)$ is surjective. This property allows one to construct an algorithm using some discretization/partition of $W$, which can be utilized to obtain an approximation of the Pareto front.

Theorem 3.1

([7, Theorem 3.1]) A point $x\in X$ is a weak efficient solution of Problem $(P)$, if and only if there exists some $w\in W$ such that $x\in Sol(P_{w}^{k})$, for all $k\in\{1,\ldots,\ell\}$. If $x\in WE(P)$, the required $w$ is given by

[TABLE]

Theorem 3.1 implies that, for a fixed $w\in W$, if

[TABLE]

then $x$ is a weak efficient point. On the other hand, if solutions of two subproblems are different from one another then a comparison is made between these solutions to eliminate a dominated point. If a new solution of a subproblem is non-dominated compared with the solutions of the other subproblems, then that solution is also weak efficient. Proposition 3.1 below states this fact more precisely and plays an important role in the weighted-constraint scalarization approach. The idea in the proposition is implemented in Algorithms 5–7 in the current article for the removal of dominated points.

Proposition 3.1

([7, Proposition 3.3])  Assume that $\exists\,w\in W$ such that $Sol(P_{w}^{j})\neq\emptyset$,$\forall j=1,\ldots,\ell$. Suppose that, for some $k\in\{1,\ldots,\ell\}$, $\exists\,\bar{x}_{k}\in Sol(P_{w}^{k})$ such that $\forall r\in\{1,2,\ldots,\ell\}$, $r\neq k$, $\exists\,\bar{x}_{r}\in Sol(P_{w}^{r})$, which satisfies

[TABLE]

Then $\bar{x}_{k}\in WE(P)$.

The feasible-value-constraint scalarization.  Introduced in [8] by the authors of the current paper, this scalarization technique has numerically been illustrated to efficiently approximate the boundary and the interior of the Pareto front even when the front is disconnected. It uses the specific expression of the weights given in Theorem 3.1, evaluated at a feasible point. Namely that with some $\hat{x}\in X$, one sets

[TABLE]

For $w\in W$, we obtain

[TABLE]

We will usually assume that a weight satisfies (4) (and hence (5) holds). The associated scalar problem is defined as

[TABLE]

For fixed $w\in W$, $\hat{x}\in X$ and $k\in\{1,\ldots,\ell\}$, define the solution set

[TABLE]

If $\bar{x}\in Sol(P_{\hat{x}}^{k})$ for some $k$, then $\bar{x}\in WE(P)$ [8, Theorem 4.6]. Note that this method is reminiscent of the $k$th-objective $\epsilon$-constraint approach, but uses a very specific choice of the vector $\epsilon$, the components of which are determined by the feasible objective function values. The relationship between these two approaches is made precise in the next section.

3.1 The relationship between $(P^{k}(\epsilon))$ and $(P_{\hat{x}}^{k})$

Given $\hat{x}\in X$, and $w$ as in (4), we show next that if $\epsilon\leq\hat{\epsilon}:=f(\hat{x})$, then the set of solutions of Problem ($P^{k}(\epsilon)$) contains the set of solutions of Problem ($P^{k}_{\hat{x}}$). Recall that (1) yields $f(x)\geq 0$ for every $x\in X$, so $\hat{\epsilon}\in\mathbb{R}^{\ell}_{+}$.

Theorem 3.2

Fix $\hat{x}\in X$, let $w$ be as in (4), and set $\hat{\epsilon}:=f(\hat{x})$. Then $Sol(P_{\hat{x}}^{k})={Sol}(P^{k}(\hat{\epsilon}))$.

Proof. Note that the objective functions of both problems are equivalent because $w\in W$. Regarding the constraints, we can use (5) and the definition of $\hat{\epsilon}$ to write $w_{k}f_{k}(\hat{x})=w_{i}f_{i}(\hat{x})=w_{i}\hat{\epsilon}_{i}$ for every $i\neq k$. Using again the fact that $w\in W$, the constraints in Problem ($P^{k}(\hat{\epsilon})$) can be equivalently written as

[TABLE]

which are the constraints of Problem ($P_{\hat{x}}^{k}$). Hence, with this choice of $\hat{\epsilon}$, Problems ($P^{k}(\hat{\epsilon})$) and ($P_{\hat{x}}^{k}$) are equivalent. In particular, they have the same solutions. Namely that ${Sol}(P^{k}(\hat{\epsilon}))=Sol(P_{\hat{x}}^{k})$ . $\Box$

3.2 The relationship between $(P_{w})$ and $(P_{\hat{x}}^{k})$

In Theorem 3.3 below, we show that those Pareto points which can be generated by the weighted-sum scalarization can also be generated by the feasible-value-constraint scalarization. The converse of this statement, however, is not true: those Pareto points which can be generated by the feasible-value-constraint approach may not be attainable by the weighted-sum method. We recall that every Pareto point is a solution of Problem $(P_{\hat{x}}^{k})$ by Theorem 3.1, but Problem $(P_{w})$ cannot generate some of the Pareto points located in a part/region of the front whose epigraph is a nonconvex set.

Theorem 3.3

Fix $w\in W$ and $\bar{x}\in Sol(P_{w})$. Then $\bar{x}\in Sol(P_{\bar{x}}^{j})$ holds for all $j$.

Proof. Let $\alpha\in W$ be as in (4)-(5) for $\hat{x}=\bar{x}$. Since $\bar{x}\in Sol(P_{w})$, we have that

[TABLE]

Let us note first that $\bar{x}$ is feasible for $(P^{j}_{\bar{x}})$ because the constraint set is defined using the weight vector $\alpha$. Assume that, for some $j\in\{1,\ldots,\ell\}$, $\bar{x}\notin Sol(P_{\bar{x}}^{j})$ . Then there exists $x\in X$ such that

[TABLE]

and

[TABLE]

where we have used the definition of $\alpha$ in the equality. From (7) and (8), we have that

[TABLE]

and

[TABLE]

where $x$ and $j$ are as in (7).

Since $w_{i}>0$, and using (9) and (10), we conclude that

[TABLE]

which contradicts (6). Therefore $\bar{x}\in Sol(P_{\bar{x}}^{j})$ for all $j$. $\Box$

Remark 3.1

If we assume that $w\geq 0$, instead of $w>0$, a reasoning similar to the one in the proof above cannot be used directly in getting the conclusion of Theorem 3.3. However, if we assume further that $Sol(P_{w})$ is a singleton for some weight $w\geq 0$, a trivial modification of the proof will work. Since the modified proof is straightforward, we do not elaborate it here.

3.3 The relationship between $(P^{k}(\epsilon))$ and $(P_{w}^{k})$

Theorem 3.1 states that every weak efficient point is the solution of the $k$th-objective weighted-constraint scalarization problem $(P_{w}^{k})$, for a suitable choice of weight. On the other hand, some weak efficient points may not be solutions of the $k$th-objective $\epsilon$-constraint problem $(P^{k}(\epsilon))$, if $\epsilon$ is too large.

Theorem 3.4

Assume that $\bar{x},\hat{x}\in WE(P)$ and

[TABLE]

If $\epsilon\geq\hat{\epsilon}:=f(\hat{x})$, then $\bar{x}\notin Sol(P^{k}(\epsilon)).$

Proof. Assume that there exists $\tilde{\epsilon}\geq\hat{\epsilon}$ such that $\bar{x}\in Sol(P^{k}(\tilde{\epsilon}))$. Thus, we can write

[TABLE]

Note that $f_{i}(\hat{x})=\hat{\epsilon}_{i}\leq\tilde{\epsilon}_{i}$ for all $i\neq k$. Hence, $x=\hat{x}$ satisfies the constraints of problem $P^{k}(\tilde{\epsilon})$. Therefore, using $x=\hat{x}$ in (12), gives

[TABLE]

contradicting (11). Therefore, for all $\epsilon\geq\hat{\epsilon}$, we must have $\bar{x}\notin Sol(P^{k}(\epsilon))$. $\Box$

Remark 3.2

Recall that Theorem 3.2 implies, with the choice of $\epsilon=f_{k}(\hat{x})$ where $\hat{x}\in X$, that the solution sets of $(P^{k}(\epsilon))$ and $(P_{\hat{x}}^{k})$ are the same. On the other hand, Theorem 3.4 asserts that, if the Pareto front has points which are weak efficient (but not just efficient) then, unless $\epsilon$ is chosen in an informed manner, the solution set of $(P^{k}(\epsilon))$ may not contain all of these weak efficient points. However, by Theorem 3.1 [7, Theorem 3.1]), the solution set of $(P_{w}^{k})$ contains the whole weak efficient set. This justifies the employment of the weighted-constraint scalarization $(P_{w}^{k})$ in the algorithms in the next section, rather than the scalarizations $(P^{k}(\epsilon))$, $(P_{\hat{x}}^{k})$, or indeed the weighted sum scalarization $(P_{w})$.

4 Algorithms

4.1 Grid generation techniques

In our study, we incorporate two types of grid generation techniques in designing the algorithms for approximating Pareto points. For details of these grid generation processes, the reader is referred to [15, 34], as well as [8, Sections 1.2 and 5].

Convex Hull of the Individual Minima (CHIM).  The CHIM grid generation technique was proposed by Das and Dennis, and employed as part of their Normal Boundary Intersection (NBI) method in [15]. The NBI method is arguably the most popular approach to constructing an approximation of the Pareto front. First, the vector of individual minima is obtained, and then the convex hull of these individual minima is generated. After a uniform discretization of the convex hull, and assigning a weight to each node of the discretization, the CHIM grid is created.

In the NBI method, given the CHIM grid along with the corresponding array of weights, a modified Pascoletti–Serafini scalarization of the problem is solved for each weight in the array to generate an approximation of the Pareto point. Geometric illustrations of the CHIM grid generation process for two and three objectives can be found in [8, Figure 1].

Sequential Boundary Generation (SBG).  The SBG grid generation technique was proposed by Mueller-Gritschneder et al. in [34]. We observe in our experiments that the SBG method works well when the CHIM grid is unable to generate the true boundary of the Pareto front. The SBG method constructs the Pareto front of the problem sequentially. In doing this, the SBG method solves linear programming problems. As a result, it requires more computational time than the CHIM method. A detailed description of the SBG technique can be found in [8, Sections 1.2 and 5].

As exemplified via test problems in [8], when the boundary of the Pareto front happens to be complicated (which is a common occasion with the three- and four-objective problems), the SBG grid is more suitable than the CHIM grid in getting a correct depiction of the boundary.

4.2 Descriptions of the algorithms

Table 1 provides a list of the algorithms we implement to solve two-, three- and four-objective problems. In the table, we list the number of objectives that an algorithm can handle, as well as the grid generation and the scalarization techniques that are utilized in that algorithm. Algorithms 1 and 2 can solve two-objective problems by using the CHIM grid. While Algorithms 3-6 can all solve three-objective problems, Algorithms 3 and 4 employ the CHIM grid generation technique, and Algorithms 5 and 6 employ the SBG grid generation technique. Algorithm 7 can handle four objectives and utilizes the SBG grid generation technique. As for the scalarization techniques, while the odd-numbered algorithms use the weighted-constraint scalarization, as in Problem (P${}_{w}^{k}$), the even numbered ones use the Pascoletti–Serafini scalarization, as in Problem (PS).

As indicated in the last column of Table 1, Algorithms 1 and 2 directly come from [7], and Algorithms 3–6 from [8]. It should be noted however that Algorithm 3 in the current paper was earlier provided in [8] as Algorithm 2 with steps written in an abridged form. Here we provide a fully expanded and explicit version of that algorithm as Algorithm 3 in Appendix A.1. Algorithm 7, on the other hand, does not appear elsewhere, and we provide its full description in Appendix A.2. We will use Algorithm 7 particularly for solving the challenging rocket injector design problem in Section 5.5.

5 Numerical Experiments

In this section, we test and compare Algorithms 1–7 that are described in Section 4, by means of two-, three- and four-objective integer and mixed-integer programming problems. The test problems 1–3 are integer programming problems. These problems are designed in such a way that the number of points in the Pareto front is finite and that they can be interpreted (or visualized easily) geometrically, so that we know the weak Pareto points prior to computations. These known weak Pareto, or weak efficient, points are referred to as being exact. One of our aims is to understand the capabilities of Algorithms 1–6 in approximating the set of exact weak Pareto points. The task of approximating a Pareto front is particularly challenging for the three- and four-objective cases. We test Algorithm 7 on a challenging real-life problem, namely the rocket injector design, which has four objective functions to minimize simultaneously.

We have coded all of Algorithms 1–7 in MATLAB, and utilized the optimization software BARON [47, 53] or SCIP [20, 1] or Bonmin [6] or IPOPT [56], for solving the single-objective integer and mixed-integer subproblems (P${}_{w}^{k}$), including the scalarized problems, in each algorithm. In particular, we used BARON in solving all the test problems. Additionally, we used SCIP, Bonmin and IPOPT for the last test problem, rocket injector design, in order to reconfirm the numerical results we obtained. All solver tolerances are set at $10^{-10}$. In the generation of the SBG grid, MATLAB’s linprog was used, with default options, for solving the associated linear programming (LP) problems. The computations have been performed on an HP ENVY 14 laptop with 4 GB RAM and core i5 at 1.6 GHz.

5.1 Test problem 1: Two objective functions

We consider, mainly for illustration purposes, a small problem with two objective functions and two integer variables:

[TABLE]

The problem above is easy to interpret: The feasible points of the above problem are those points with integer coordinates in the lower-left quarter of a circle of radius 4 centred at $(4,4)$. One can easily count that there are 17 feasible points of the problem. Then, of these feasible points, the ones which are located as far to the “south-west” as possible are weak Pareto points. There are nine weak Pareto points, the set of which can simply be written as (see Figure 1(a))

[TABLE]

The above set is nothing but the Pareto front, which is discrete and so is disconnected. Note that $f_{1}(x)=x_{1}$ and $f_{2}(x)=x_{2}$ are the objective function of the problem. As the utopia point, we take $u=(-10,-10)$. For Algorithms 1 and 2, we use the CHIM grid. Since the CPU time requirements are in general different for each of the algorithms for the same number of the grid points, we have adjusted the number of grid points in such a way that the CPU time each algorithm takes is roughly the same (in this case, around 13 and 14 seconds, respectively). This way, we provide the two algorithms an equal footing in the comparisons. The data on the CPU time, the number of subproblems attempted, and the number of weak Pareto points generated by each algorithm, are reported in Table 2. Parts (a) and (b) of Figure 1 display all of the Pareto points, i.e., an approximation of the Pareto front, that could be obtained by each algorithm, respectively, within the CPU times indicated.

Algorithm 1, as shown by Figure 1(a), successfully finds all nine weak Pareto points in the Pareto front. On the other hand, Algorithm 2 is able to generate only 6 of these points. When we increase the number of grid points, i.e., increase the allowed CPU time, by (say) ten times, the number of Pareto points obtained by Algorithm 2 still turns out to be 6. In other words, no improvement in the approximation of the Pareto front can be obtained by Algorithm 2, even after providing the algorithm with a much finer grid and so allowing a much longer computational time.

5.2 Test problem 2: Three objective functions – nonlinear constraints

Next we modify the test problem 1 in Section 5.1 by introducing one more objective function and one more variable as follows.

[TABLE]

The feasible region is now defined to be a spherical ball rather than a part of a circular region; namely, the feasible set consists of the points with integer coordinates in a sphere of radius 2 centred at $(2,2,2)$. In this case, the problem has 70 feasible points, of which 19 are weak Pareto. The set of weak Pareto points, or the Pareto front, can simply be written as (see Figures 2(a) and 2(c))

[TABLE]

We note that $f_{1}(x)=x_{1}$, $f_{2}(x)=x_{2}$ and $f_{3}(x)=x_{3}$. We choose the utopia point as $u=(-10,-10,-10)$. In order to solve the problem, i.e., to obtain an approximation of the Pareto front, we have implemented Algorithms 3 and 4, both of which use the CHIM grid. As in the case of Test Problem 1, with the choice of the same CHIM grid, each algorithm requires a substantially different length of CPU time. Therefore, we have adjusted the number of points in the CHIM grids of the algorithms in such a way that the CPU time each algorithm takes is roughly the same (in this case, about 150 seconds each), so that the two algorithms can be compared on an equal footing. Table 3 summarizes the performance of each algorithm.

As can be seen in Table 3, Algorithm 4 is able to find only 7 weak Pareto points out of 19, as can also be seen in Figures 2(b) and (d). Note that the Pareto front in Figure 2(d) is a rotated view of the same front in Figure 2(b) and has been added to the figure to make the visualization easier. Next we provided Algorithm 4 with much finer CHIM grids, which obviously resulted in far longer computational times; however, not a single new weak Pareto point could be obtained by Algorithm 4. This indicates that Algorithm 4 is not efficient for this test problem. On the other hand, Algorithm 3 generates all the 19 weak Pareto points (see Figures 2(a) and 2(c)) in 147 seconds of CPU time.

5.3 Test problem 3: Three objective functions – linear constraints

This time, we consider the test problem given by Antunes et al. [2, Exercise 5, Section 6.5], who studied certain properties of two known weak Pareto points, but not the whole Pareto front. We are interested here in constructing the whole Pareto front of this problem, which, to the best of our knowledge, has not been done elsewhere, yet.

[TABLE]

The problem has three non-negative integer decision variables. The feasible region is defined by the intersection of four closed half-spaces as indicated above and has 83 feasible points (with non-negative integer coordinates), of which 60 are weak Pareto. The (disconnected) Pareto front consists of these 60 discrete points – see the graphs in parts (a) and (c) of Figure 4.

The three objective functions of the problem are namely $f_{1}(x)=-x_{1}$, $f_{2}(x)=-x_{2}$ and $f_{3}(x)=-x_{3}$. We choose $u=(-100,-100,-100)$ as the utopia point. For solving this problem, we used Algorithms 3–6. For the features of these algorithms, see Table 1. As in the test problems 1 and 2, for fairness in comparisons, we adjusted the number of points in, depending on the algorithm, the CHIM or SBG grids, in such a way that the CPU time each algorithm takes is roughly the same (around 210 seconds for Algorithms 3–4 and 130 seconds for Algorithms 5–6). Table 4 summarizes the numerical performance of Algorithms 3 and 4, while Figures 3 and 4 show the Pareto points generated by each algorithm.

For the given CHIM grid, or the allowed CPU time of 202 seconds, Algorithm 3 generates 58 of the 60 weak Pareto points in the front. On the other hand, Algorithm 4 is able to find only around half of the weak Pareto points, in about the same CPU time as that of Algorithm 3. The Pareto fronts generated by these algorithms are depicted in Figure 3. When we provide a much finer CHIM grid, i.e., increase the allowed CPU time by about 10 fold, either algorithm fails to generate any new points in the front. We conclude that Algorithm 3 is more efficient and powerful than Algorithm 4 in approximating the Pareto front.

Algorithm 5, which implements the SBG grid, generates all of the 60 weak Pareto points in the front in just 128 seconds, whereas Algorithm 6 produces only 33 out of the 60 Pareto points, in about the same amount of time. When we provide a much finer SBG grid, resulting in a CPU time of 120 minutes or longer, Algorithm 6 can find at most 44 of the weak Pareto points. Considering all of Algorithms 3–6, we conclude that Algorithm 5 is the most efficient and powerful.

To test performance when some of the variables are continuous, we transform Test Problem 3 into a mixed-integer problem in two separate instances:

1. (i)

$x_{1}$ and $x_{2}$ are continuous, and $x_{3}$ an integer. 2. (ii)

$x_{2}$ is continuous, and $x_{1}$ and $x_{3}$ integers.

We choose $u=(-100,-100,-100)$ as the utopia point. To solve either of the instances (i) and (ii), we use both of Algorithms 3 and 4 and compare. In order to make comparisons on an equal footing, we adjust the number of CHIM points in the algorithms in such a way that the CPU time each algorithm takes is roughly the same. For case (i), we set around 500 seconds, and for case (ii), around 300 seconds.

The Pareto fronts obtained for Instance (i) are shown in Figures 5(a) and (b) for Algorithms 3 and 4, respectively. We observe that while Algorithm 4 cannot approximate the upper part of the Pareto front as can be seen in Figure 5(b), Algorithm 3 efficiently approximates all six plane segments of the Pareto front as shown in Figure 5(a).

In Instance (ii), we report a similar kind of performance: As can be observed in Figures 5(c) and (d), respectively, Algorithm 3 constructs a better approximation of the Pareto front than that by Algorithm 4. Algorithm 4, which employs the Pascoletti–Serafini scalarization, clearly misses many parts of the front in either instance.

5.4 Test problem 4: Three objective functions – nonlinear constraints

We consider next a test problem given by de Santis et al. [16, Test Instance 4.5, Section 4.3], which involves three objective functions and integer, as well as continuous, variables, as stated below.

[TABLE]

The challenge here is to construct (approximately) the whole Pareto front of this problem, and observe the performances of Algorithms 3–6 in carrying out this task. As can be seen, the problem has three continuous variables and one integer variable. The resulting Pareto front consists of five disconnected Pareto surfaces—see Figure 6. We chose $u=(-100,-100,-100)$ as the utopia point. We adjusted the number of SBG grid points in the algorithms in such a way that the CPU times needed by the algorithms were more or less the same. Therefore, their performances could be compared on a more equal footing. For Algorithms 5 and 6, we set around 1200 seconds of CPU time. Both algorithms generated identical Pareto fronts, which are depicted in Figure 6.

We also implemented Algorithms 3 and 4 to solve Test Problem 4. Although Algorithm 3 can be seen to yield a better approximation of the Pareto front than Algorithm 4 (see Figure 7), this approximation is still nowhere nearly as good as the one obtained by Algorithms 5–6. It is conceivable to think that the reason Algorithms 5 and 6 do a better job is that they employ the SBG grid, which is particularly effective when the Pareto front has a complicated boundary, as is the case in this example problem.

Remark 5.1

We aim to generate an approximation of the entire Pareto front, i.e., including also all of the weak efficient points, and we stress that aiming only for the efficient points may miss completely some sizeable and significant portions of the Pareto front. We can elaborate on the significance and justification of this view point as follows.

In the mixed-integer version of Test Problem 3, we find that for any (fixed) assignment of $x_{3}$ and $x_{1}$ there is exactly one choice/value of $x_{2}$ that leads to an efficient solution. By looking at Figure 5(c), if we aim only for the efficient points but not the weak efficient ones, any of the line segments appearing in the (weak) Pareto front would be reduced to a single point, leading to a grossly inaccurate approximation of the front.

While our algorithms are successful in obtaining weak efficient solutions which are not efficient, they are also successful in obtaining efficient solutions, when compared with other algorithms. For example, while Algorithm 3 can find all the efficient points in both instances (i) and (ii) of Problem 3—see Figure 5(a) and (c)—Algorithm 4 fails in doing so.

5.5 An application to rocket injector design

The liquid-rocket injector design problem was previously studied as a multi-objective optimization problem in [8, 21, 55]. Two primary objectives are of concern in this design problem: (i) improvement of the performance of the injector and (ii) increase of its survivability or lifetime. The performance of the injector is influenced by the axial length of the thrust chamber while the lifetime of the injector is associated with the thermal field inside the thrust chamber. For an illustration and visual representation of the injector design we refer the reader to [21]. There is a conflicting interplay between these two main objectives: while high temperatures improve the performance, they reduce the lifetime and thus decrease survivability. Four design variables are introduced in [21] to construct the mathematical model of the rocket injector design problem; namely,

• $x_{1}$:

hydrogen flow angle,

• $x_{2}$:

oxidizer post tip thickness,

• $x_{3}$:

decrement with respect to the baseline cross-section area of the tube carrying oxygen,

• $x_{4}$:

increment with respect to the baseline cross-section area of the tube carrying hydrogen.

In our present study, we introduce an integer-valued variable $\tilde{x}_{1}\in\{0,1,2,3\}$ and impose a new constraint $x_{1}=0.2\,\tilde{x}_{1}$. The constraint enforces the hydrogen flow angle ($x_{1}$) to be one of the prescribed set of angles $\{0,0.2,0.4,0.6\}$. We note that, with the introduction of an integer variable, the modified problem has a disconnected domain and so is expected to have a disconnected Pareto front.

We consider four objective functions as in [21]; namely,

$f_{1}(x)$:

face temperature, which is the maximum temperature of the injector face,

$f_{2}(x)$:

tip temperature, which is the maximum temperature on the post tip of the injector,

$f_{3}(x)$:

combustion length, which is the distance from the inlet where 99% of the combustion are complete,

$f_{4}(x)$:

wall temperature, which is the wall temperature at three inches (fourth probe) from the injector face.

The modified rocket injector design problem considered in this paper is a multi-objective mixed-integer optimization problem and described as follows.

[TABLE]

where

[TABLE]

subject to

[TABLE]

[TABLE]

In obtaining an approximation of the Pareto front of the above problem in the four-dimensional value space, we have used the new Algorithm 7, which implements the weighted-constraint scalarization and the SBG grid. Note that an implementation of the CHIM grid generation technique (instead of the SBG grid) in Algorithm 7 would not have worked for the rocket injector design problem, because of the rather complex boundary of the Pareto front of this problem. Hence Algorithm 7 uses the SBG grid.

We compare our mixed-integer multi-objective programming results with those obtained for the case of continuous variables in [8]. We choose to display the projection $(f_{1},f_{2},f_{3},f_{4})\mapsto(f_{1},f_{2},f_{3},0)$ in the $f_{1}f_{3}f_{2}$-space for this comparison. The Pareto front depicted in parts (a) and (c) of Figure 8 corresponds to the case of continuous variables in [8], and was obtained by using Algorithm 9 from [8]. On the other hand, the Pareto front in parts (b) and (d) of Figure 8 corresponds to the case of mixed integer variables, and have been obtained by using Algorithm 7.

Algorithm 7 obtained $41014$ points, which, after the elimination of dominated points, resulted in $33536$ points representing the Pareto front. Note that what we obtain is actually a discretization of the front, i.e., the continuous parts of the front are approximated by, in particular computed as, a discrete set of points. A computed point may in fact not be a weak Pareto point due to computational errors. Therefore, it is possible to classify such erroneous points as non-dominated when compared with the other computed points. For this reason, one cannot ensure that all dominated points have been eliminated, or weeded out; however, this example application still illustrates the utility of our approach. The Pareto front approximation in Figure 8(a) for continuous variables had taken about 2 hours, using the same computer, as reported in [8]. The Pareto approximation in Figure 8(b), on the other hand, took about 32 hours. One should recall that solving mixed integer problems in general requires substantially more computer memory and CPU time.

5.5.1 Comparison of Algorithm 7 with a brute force or slicing technique

In order to obtain an approximation of the Pareto front, we aim to propose algorithms supported with mathematical results rather than a brute-force methodology. The following reasons motivate our approach for not resorting to brute-force.

• •

The computational time required by a brute-force algorithm will depend exponentially on the number of optimization variables and the number of objective functions.

• •

A more accurate representation of the front would typically require an exponentially growing number of points and effort, which might grind the process down to a halt.

• •

If the Pareto front’s boundary is fragmented or disconnected, i.e., if the front itself has gaps on the “inside of itself,” then one would need even more points to approximate the Pareto front. Many of the points generated by algorithms around these gaps are typically non-Pareto points and, therefore, require more computational time. This case creates a bigger concern for the brute-force algorithm.

Here we aim to justify our choice of not resorting to brute force, or the so-called slicing technique, by carrying out experiments on the rocket injector design problem with each of the methods (brute and non-brute) and making comparisons. Both methods are allowed to take about (the same) 2.8 hours of CPU time. Here we first run Algorithm 7 for 2.8 hours with slicing, i.e., run Algorithm 7 four times consecutively, each time with a fixed value of $x_{1}=0,0.2,0.4,0.6$, and superpose all the solutions and do the necessary weeding (or filtering) out. We run Algorithm 7 again for 2.8 hours, but this time only once, without fixing the value of $x_{1}$. Ultimately we compare the performances of each approach by looking the Pareto fronts that they generate and see which of these Pareto fronts come closer to the refined/accurate front presented in Figures 8(b) and 8(d).

In the approximation of the Pareto front depicted in Figure 8(b) and 8(d), we had used the mixed-integer problem solver SCIP for the subproblems of Algorithm 7. On the other hand, when we apply the slicing technique we get a continuous optimization problem for every slice, i.e., for every fixed value of $x_{1}$. Therefore we use a continuous solver, namely Ipopt with a MATLAB interface, for the slicing technique, in order to exploit fully the computational power of a continuous solver. It is well-known that continuous solvers are far more efficient than mixed-integer solvers in solving a continuous problem, in general.

In running Algorithm 7 for the rocket injector design problem, we have adjusted the number of SBG grid points in Algorithm 7 so as to adjust the CPU time to be about $2.8$ hours. For each of the four fixed values of $x_{1}$, Algorithm 7 has generated four Pareto fronts for the four sliced-up problems. Adding all the points in these fronts has then yielded $21611$ “candidate” solutions. Discarding/weeding out the dominated points leaves over a mere $8571$ non-dominated points, which are shown in Figures 9(a) and 9(c). This illustrates that only about 40% of the points found by Algorithm 7 via slicing are non-dominated.

We have also solved the mixed-integer rocket injector design problem by running Algorithm 7 without any slicing, i.e., without fixing the value of $x_{1}$. In this case, we have set the same 2.8 hours of CPU time by providing the size of the SBG grid to Algorithm 7 accordingly. Since the subproblems of Algorithm 7 (without slicing) are mixed-integer problems, we have tried each of the mixed-integer solvers Bonmin, BARON and SCIP, and observed that Bonmin was the fastest. With the choice of Bonmin for solving the subproblems, Algorithm 7 generated $21882$ candidate points. Of these points, after weeding out, we obtained $19698$ non-dominated points, which are about 90% of all points computed. In this particular example, the number of Pareto points found by Algorithm 7 without slicing is about 2.3 times higher than that with slicing.

The Pareto front obtained by Algorithm 7 without slicing is depicted in Figures 9(b) and 9(d) from two distinct viewing angles. These viewing angles are the same as those used in Figures 8(b) and 8(d), and in 9(a) and 9(c), respectively, for ease of comparisons. First, we observe that a triangle-shaped part of the front is missing from the approximation of the front presented in Figure 9(c)—compare with Figures 8(d) and 9(d). We also observe that the boundary of the front in Figures 9(b) and 9(d) look better developed and defined. Overall, the images of the front in Figures 9(b) and 9(d) seem to be closer in appearance to those in Figures 8(b) and 8(d), compared with those in Figures 9(a) and 9(c). Therefore, based on the numerical experiments and results for this particular example, we conclude that Algorithm 7 without slicing does a better job than that with slicing.

We note that filtering, or weeding out, was also carried out for all the test examples other than the rocket injector design problem, but no dominated points were found in those cases. Filtering seems to be more important, or necessary, for the rocket injector design problem, since this problem is much more difficult than the previous ones. There are cases in the rocket injector design problem when a point computed by the optimization software is only locally Pareto optimum. In such a case, we need to weed out these points, if we can, by making comparisons with the other points computed. However, if one has not found many other Pareto points around a “local” (dominated) Pareto point, then it may not be possible to remove that dominated only-locally-Pareto point. We have to accept, on the other hand, that this situation is an inherent deficiency/property of any general algorithm, unless the problem in question has a particular structure that can be exploited to find global Pareto optimal points.

6 Conclusion and Discussion

We have established relationships amongst various scalarization techniques listed in Section 3. We have implemented and tested algorithms, which were originally proposed for continuous multi-objective problems, for constructing the Pareto front of integer and mixed-integer programming problems. These algorithms utilize one of the CHIM and SBG grid generation techniques and one of the Pascoletti–Serafini or weighted-constraint scalarization techniques, with two or three objective functions to minimize. We have also presented a new four-objective algorithm. We applied the algorithms to test problems, one of them being a challenging four objective mixed-integer programming problem.

We found that for the given test problems, (i) Algorithm 1 outperforms Algorithm 2 for the two-objective test problem, (ii) Algorithm 3 outperforms Algorithm 4 for the first three-objective test problem and (iii) Algorithm 5 outperforms Algorithm 6 for the second three-objective test problem. The common characteristic of the algorithms which are successful is that they use the weighted-constraint scalarization. This scalarization approach was previously demonstrated in [7, 8] to perform well for continuous-variable problems with disconnected feasible domains. It is also worth noting that, in solving the second three-objective test problem, we have demonstrated that Algorithm 5, which utilizes the SBG grid and the weighted-constraint scalarization, outperforms not only Algorithm 6 but also Algorithms 3 and 4.

We have used the weighted-constraint scalarization and the SBG grid in Algorithm 7, which we have presented in this paper for the first time – see Appendix A.2. We have proposed a modified, mixed-integer programming, version of the already challenging rocket injector design problem (previously studied as a continuous problem in [8]), and applied Algorithm 7 to construct the Pareto front of this problem successfully. The choice of the SBG grid for Algorithm 7 is essential as the boundary of the Pareto front of this problem is quite complicated.

The continuous version of the rocket injector design problem was solved in [8] with an algorithm (Algorithm 9 in [8]) implementing the feasible-value-constraint scalarization (also introduced in [8]) and the SBG grid. In the same paper, an algorithm with the same scalarization and the CHIM grid (instead of the SBG grid) was not successful. In view of Remark 3.2 and the particular success of the weighted-constraint scalarization with problems with disconnected domains, we have proposed in the current paper Algorithm 7 implementing this scalarization along with the SBG grid.

Acknowledgments

The authors would like to thank Jörg Fliege once again for pointing to the challenging multi-objective rocket injector design problem, which was earlier used by the authors in [8]. In the present paper, the rocket injector problem has been modified/altered so as to make it a multi-objective mixed-integer problem, yielding an even more challenging instance.

Appendix A

A.1  Algorithm 3

We use the CHIM grid and the weighted-constraint scalarization in the algorithm. In Step 2 of Algorithm 3, each objective function, subject to the original constraints of the problem, is minimized. These individual minima are used in Step 3 to form a triangle shaped grid. Each grid point corresponds to a weight vector in $\mathbb{R}^{3}$. In Step 4, three sub-problems are solved at each grid point to generate Pareto points.

In Step 4(b), the approximation of the weak efficient point is calculated using the fact that if $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}=:\bar{x}\;\mbox{(say)}$ holds, then the solution $\bar{x}$ would be weak efficient. Here, $\bar{x}_{k}$ are the solutions of (P${}_{w}^{k}$) for $k=1,...,3$. On the other hand, if $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}$ does not hold, then any dominated point is removed from the set {$\bar{x}_{1}$, $\bar{x}_{2}$, $\bar{x}_{3}$} (see Step 4(b(ii))). The latter case is typically encountered when the Pareto front and/or the domain is disconnected. Therefore, this algorithm is efficient in finding Pareto points even when the feasible set is discrete or disconnected.

Step $\mathbf{1}$

(Input)

Choose the utopia point $u=(u_{1},u_{2},u_{3})$ and $x\in X$. Set the number of partition points $N$ in the $(N+1)(N+2)/2$ grid points. Set $s=0$.

Step $\mathbf{2}$

**(Determine the individual minima)

**Solve Problem $(P_{i}):=\displaystyle\min_{x\in X}f_{i}(x),\;\;i=1,2,3$, that give the solutions $\bar{x}_{f_{i}}$, for $i=1,2,3$, respectively.

Set $\bar{F}:=[f_{1}(\bar{x}_{f_{3}}),f_{2}(\bar{x}_{f_{3}}),f_{3}(\bar{x}_{f_{3}})]$

Step $\mathbf{3}$

(Generate grid points over the CHIM grid)

Compute $a_{k}$, $b_{k}$ and $c_{k}$, for $k=1,2$, using [8, Sections 5.1–5.2] : Let

[TABLE]

for $i=1,2,3$. Then

[TABLE]

For $i_{1}:=0,1,...,(N-1)$

{

        Let $\hat{w}_{k}^{i_{1}}:=a_{k}+i_{1}(c_{k}-a_{k})/N,\;\;k=1,2,$

        and $\tilde{w}_{r}^{i_{1}}:=b_{r}+i_{1}(c_{r}-b_{r})/N,\;\;r=1,2,$

        For $i_{2}:=0,1,...,(N-i_{1})$

        {

              Let $w_{\ell}^{(i_{1},i_{2})}:=\hat{w}_{\ell}^{i_{1}}+i_{2}(\tilde{w}_{\ell}^{i_{1}}-\hat{w}_{\ell}^{i_{1}})/(N-i_{1}),\;\;\ell=1,2,$

        }

}

Step $\mathbf{4}$

**(Solve scalar sub-optimization problems)

**For $i_{1}:=0,1,...,(N-1)$

{

        For $i_{2}:=0,1,...,(N-i_{1})$

        {

              Let $w_{k}=w_{k}^{(i_{1},i_{2})},\;\;k=1,2.$ Let $w=(w_{1},w_{2},1-w_{1}-w_{2})$.

(a)

Find $x_{k}$ that solves Problems ($P_{w}^{k}$), $k=1,2,3.$

(b)

Determine weak efficient points :

Let $s:=s+1.$

(i)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}$, then set $\bar{x}=\bar{x}_{1}$ (a weak efficient point)

and $F(s):=[f_{1}(\bar{x}),f_{2}(\bar{x}),f_{3}(\bar{x})].$

(ii)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}$ does not hold, then, any dominated point is discarded by the following Steps 1-3.

(1)

If $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{1})$ and $f_{3}(\bar{x}_{3})\geq f_{3}(\bar{x}_{1})$ then, $\bar{x}_{1}$ weak efficient point, and $F(s):=[f_{1}(\bar{x}_{1}),f_{2}(\bar{x}_{1}),f_{3}(\bar{x}_{1})].$

(2)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{2})$ and $f_{3}(\bar{x}_{3})\geq f_{3}(\bar{x}_{2})$ then, $\bar{x}_{2}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{2}),f_{2}(\bar{x}_{2}),f_{3}(\bar{x}_{2})].$

(3)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{3})$ and $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{3})$ then, $\bar{x}_{3}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{3}),f_{2}(\bar{x}_{3}),f_{3}(\bar{x}_{3})].$

}

}

Step $\mathbf{5}$

(Output)

Set $F(s+1):=\bar{F}$.

The array of Pareto points $F$, is an approximation of the Pareto front.

A.2  Algorithm 7

This algorithm is an extension of Algorithm 5 to the four-objective case. In Algorithm 7, the Pareto front is approximated in six main steps. In Step 2, the individual minima of all four objectives are obtained. In Step 3, the boundaries of the Pareto front of the two-objective subproblems are constructed. In Step 4, a linear programming problem is solved, to obtain base points that are used in the construction of the SBG grid in the spaces of the objective function triplets. In Step 5, the boundaries of the Pareto fronts of the three-objective subproblems are generated. In Step 6, another linear programming problem is solved, to obtain base points that are used in the construction of the SBG grid in the space of all objective function quadruplets. Finally, in Step 7, a four-objective minimization problem is solved and so the points within the boundary of the Pareto front are generated.

Step 1

(Input)

Set the number of partition points $N$.

Step 2

(Determine the individual minima)

Find $\bar{x}_{f_{k}}$ that solves Problem (Pk):= $\displaystyle\min_{x\in X}f_{k}$, for $k=1,2,3,4$.

Set $F(k):=[f_{1}(\bar{x}_{f_{k}}),f_{2}(\bar{x}_{f_{k}}),f_{3}(\bar{x}_{f_{k}}),f_{4}(\bar{x}_{f_{k}})]$ and $W(k):=e_{k}$, for $k=1,2,3,4$.

Set $s=k+1$.

Step 3

(Generate Pareto fronts of the objective pairs)

Set the direction vector $v^{(i,k)}$ in the $f_{i}f_{k}$-space, $i=1,2,3$, $k=2,3,4$, $i\neq k$ :

Let $v^{(i,k)}:=[v_{1}^{(i,k)},v_{2}^{(i,k)}]$.

Set $v_{1}^{(i,k)}=\max\{f_{i}(\bar{x}_{f_{i}}),f_{i}(\bar{x}_{f_{k}})\}-f_{i}(\bar{x}_{f_{i}})$

and $v_{2}^{(i,k)}:=\max\{f_{k}(\bar{x}_{f_{i}}),f_{k}(\bar{x}_{f_{k}})\}-f_{k}(\bar{x}_{f_{k}})$ .

Auxiliary weights for a parallel ray in the $f_{i}f_{k}$-space :

Let $c=1/(v_{1}^{(i,k)}+v_{2}^{(i,k)})$, $\tilde{w}_{1}^{(i,k)}=cv_{1}^{(i,k)}$, $\tilde{w}_{2}^{(i,k)}=1-\tilde{w}_{1}$.

For $j:=1,...,(N-1)$

{

(a)

(Generate Pareto points for $\{f_{1},f_{2}\}$)

Let $w_{1}=(N-j)/N$; $w_{2}=1-w_{1}$; $w_{3}=0$; $w_{4}=0$;

Set $W(s):=[w_{1},w_{2},w_{3},w_{4}]$ and $s=s+1$.

Form a base point for $\{f_{1},f_{2}\}$:

Let $u_{1}=w_{1}f_{1}(\bar{x}_{f_{1}})+w_{2}f_{1}(\bar{x}_{f_{2}})$; $u_{2}=w_{1}f_{2}(\bar{x}_{f_{1}})+w_{2}f_{2}(\bar{x}_{f_{2}})$.

(i)

Find $\bar{x}_{1}$ and $\bar{x}_{2}$ that solve Problems (P${}_{w}^{1}$) and (P${}_{w}^{2}$), respectively.

(iii)

Determine weak efficient points :

• $\bullet$

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}$ (a weak efficient point), then set

$F(s):=[f_{1}(\bar{x}),f_{2}(\bar{x}),f_{3}(\bar{x}),f_{4}(\bar{x})]$ and $s=s+1$.

• $\bullet$

If $\bar{x}_{1}\neq\bar{x}_{2}$, then $\bar{x}_{1}$ and $\bar{x}_{2}$ are weak efficient points if they are not dominated by each other. The following steps remove a dominated point, if it exists.

(1)

If $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{1})$ then $\bar{x}_{1}$ weak efficient point and

$F(s):=[f_{1}(\bar{x}_{1}),f_{2}(\bar{x}_{1}),f_{3}(\bar{x}_{1}),f_{4}(\bar{x}_{1})].$

(2)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{2})$ then $\bar{x}_{2}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{2}),f_{2}(\bar{x}_{2}),f_{3}(\bar{x}_{2}),f_{4}(\bar{x}_{2})],$ and $W(s+1)=W(s)$.

(b)

(Generate Pareto points for $\{f_{2},f_{3}\}$)

Let $w_{1}=0$; $w_{2}=(N-j)/N$; $w_{3}=1-w_{2}$; $w_{4}=0$.

Do the same as in the rest of Step 3(a) for the new objective pair $\{f_{2},f_{3}\}$ instead of $\{f_{1},f_{2}\}$.

(c)

(Generate Pareto points for $\{f_{3},f_{4}\}$)

Let $w_{1}=0$; $w_{2}=0$; $w_{3}=(N-j)/N$; $w_{4}=1-w_{3}$.

Do the same as in the rest of Step 3(a) for the new objective pair $\{f_{3},f_{4}\}$ instead of $\{f_{1},f_{2}\}$.

(d)

(Generate Pareto points for $\{f_{1},f_{3}\}$)

Let $w_{1}=(N-j)/N$; $w_{2}=0$; $w_{3}=1-w_{1}$; $w_{4}=0$.

Do the same as in the rest of Step 3(a) for the new objective pair $\{f_{1},f_{3}\}$ instead of $\{f_{1},f_{2}\}$.

(e)

(Generate Pareto points for $\{f_{1},f_{4}\}$)

Let $w_{1}=(N-j)/N$; $w_{2}=0$; $w_{3}=0$; $w_{4}=1-w_{1}$.

Do the same as in the rest of Step 3(a) for the new objective pair $\{f_{1},f_{4}\}$ instead of $\{f_{1},f_{2}\}$.

(f)

(Generate Pareto points for $\{f_{2},f_{4}\}$)

Let $w_{1}=0$; $w_{2}=(N-j)/N$; $w_{3}=0$; $w_{4}=1-w_{2}$.

Do the same as in the rest of Step 3(a) for the new objective pair $\{f_{2},f_{4}\}$ instead of $\{f_{1},f_{2}\}$.

}

Step 4

(Form base points in the spaces of objective triplets)

Solve an LP problem as in [34, Section 7.4] to get an array of base points

$U^{(1,2,3)}(r)=[u_{1}(r),u_{2}(r),u_{3}(r),u_{4}(r)]$, $r=1,\ldots,(N-1)(N-2)/2$.

Form base points $U^{(2,3,4)}(r)$, $U^{(1,3,4)}(r)$ and $U^{(1,2,4)}(r)$ in the $f_{2}f_{3}f_{4}$-, $f_{1}f_{3}f_{4}$- and $f_{1}f_{2}f_{4}$-spaces, respectively.

Step 5

(Generate Pareto fronts of the objective triplets)

(a)

(Generate interior Pareto points for $\{f_{1},f_{2},f_{3}\}$)

Set the direction vector $v$ and weights for a parallel ray in the $f_{1}f_{2}f_{3}$-space:

Let $v:=[v_{1},v_{2},v_{3}]$. Set $v_{j}:=\max\{f_{j}(\bar{x}_{f_{1}}),f_{j}(\bar{x}_{f_{2}}),f_{j}(\bar{x}_{f_{3}})\}-f_{j}(\bar{x}_{f_{j}})$, $j=1,2,3$.

Let $c=1/(v_{1}+v_{2}+v_{3})$, $w_{1}=cv_{1}$, $w_{2}=cv_{2}$ and $w_{3}=1-w_{1}-w_{2}$.

For $r=1,\ldots,\mbox{length}(U^{(1,2,3)}(r))$

{

(a)

Find $x_{k}$ that solves Problems ($P_{w}^{k}$), $k=1,2,3.$

(b)

Determine weak efficient points :

Let $s:=s+1.$

(i)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}$ then set, $\bar{x}=\bar{x}_{1}$ (a weak efficient point)

and $F(s):=[f_{1}(\bar{x}),f_{2}(\bar{x}),f_{3}(\bar{x}),f_{4}(\bar{x})].$

(ii)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}$ does not hold, then, do the same as in Steps 4(b[ii]) of Algorithm 1.

}

(b)

(Generate interior Pareto points for $\{f_{2},f_{3},f_{4}\}$, $\{f_{1},f_{3},f_{4}\}$ and $\{f_{1},f_{2},f_{4}\}$)

Do the same as in the rest of Step 5(a) for the new objective triplets $\{f_{2},f_{3},f_{4}\}$, $\{f_{1},f_{3},f_{4}\}$ and $\{f_{1},f_{2},f_{4}\}$ instead of $\{f_{1},f_{2},f_{3}\}$, respectively.

Step 6

(Form base points in the $f_{1}f_{2}f_{3}f_{4}$-space)

Given the Pareto points obtained in Steps 3 and 5 solve an LP problem to construct base points in the $f_{1}f_{2}f_{3}f_{4}$-space.

Step 7

(Generate interior Pareto points for $\{f_{1},f_{2},f_{3},f_{4}\}$)

Set the direction vector $v$ and weights for a parallel ray in the $f_{1}f_{2}f_{3}f_{4}$-space:

Let $v:=[v_{1},v_{2},v_{3},v_{4}]$. Set $v_{j}:=\max\{f_{j}(\bar{x}_{f_{1}}),f_{j}(\bar{x}_{f_{2}}),f_{j}(\bar{x}_{f_{3}}),f_{j}(\bar{x}_{f_{4}})\}-f_{j}(\bar{x}_{f_{j}})$,

$j=1,2,3,4$.

Let $c=1/(v_{1}+v_{2}+v_{3}+v_{4})$, $w_{1}=cv_{1}$, $w_{2}=cv_{2}$ , $w_{3}=cv_{3}$ and $w_{4}=1-w_{1}-w_{2}-w_{3}$.

For $r=1,\ldots,\mbox{length}(U^{(1,2,3,4)}(r))$

{

(a)

Find $x_{k}$ that solves Problems ($P_{w}^{k}$), $k=1,2,3,4.$

(b)

Determine weak efficient points :

Let $s:=s+1.$

(i)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}=\bar{x}_{4}$ then set, $\bar{x}=\bar{x}_{1}$ (a weak efficient point)

and $F(s):=[f_{1}(\bar{x}),f_{2}(\bar{x}),f_{3}(\bar{x}),f_{4}(\bar{x})].$

(ii)

If $\bar{x}_{1}=\bar{x}_{2}=\bar{x}_{3}=\bar{x}_{4}$ does not hold, then, any dominated point is discarded by the following Steps 1-4.

(1)

If $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{1})$ and $f_{3}(\bar{x}_{3})\geq f_{3}(\bar{x}_{1})$ and $f_{4}(\bar{x}_{4})\geq f_{4}(\bar{x}_{1})$ then, $\bar{x}_{1}$ weak efficient point, and $F(s):=[f_{1}(\bar{x}_{1}),f_{2}(\bar{x}_{1}),f_{3}(\bar{x}_{1}),f_{4}(\bar{x}_{1})].$

(2)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{2})$ and $f_{3}(\bar{x}_{3})\geq f_{3}(\bar{x}_{2})$ and $f_{4}(\bar{x}_{4})\geq f_{4}(\bar{x}_{2})$ then, $\bar{x}_{2}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{2}),f_{2}(\bar{x}_{2}),f_{3}(\bar{x}_{2}),f_{4}(\bar{x}_{2})].$

(3)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{3})$ and $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{3})$ and $f_{4}(\bar{x}_{4})\geq f_{4}(\bar{x}_{3})$ then, $\bar{x}_{3}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{3}),f_{2}(\bar{x}_{3}),f_{3}(\bar{x}_{3}),f_{4}(\bar{x}_{3})].$

(4)

If $f_{1}(\bar{x}_{1})\geq f_{1}(\bar{x}_{4})$ and $f_{2}(\bar{x}_{2})\geq f_{2}(\bar{x}_{4})$ and $f_{3}(\bar{x}_{3})\geq f_{3}(\bar{x}_{4})$ then, $\bar{x}_{4}$ weak efficient point. Let $s:=s+1,$ and $F(s):=[f_{1}(\bar{x}_{4}),f_{2}(\bar{x}_{4}),f_{3}(\bar{x}_{4}),f_{4}(\bar{x}_{4})].$

}

Step 8

(Output)

The array of Pareto points, $F$, is an approximation of the Pareto front.