From a Pareto Front to Pareto Regions: A Novel Standpoint for Multiobjective Optimization

Abstract

This work presents a novel approach for multiobjective optimization problems, extending the concept of a Pareto front to a new idea of the Pareto region. This new concept provides all the points beyond the Pareto front, leading to the same optimal condition with statistical assurance. This region is built using a Fisher–Snedecor test over an augmented Lagragian function, for which deductions are proposed here. This test is meant to provide an approximated depiction of the feasible operation region while using meta-heuristic optimization results to extract this information. To do so, a Constrained Sliding Particle Swarm Optimizer (CSPSO) was applied to solve a series of four benchmarks and a case study. The proposed test analyzed the CSPSO results, and the novel Pareto regions were estimated. Over this Pareto region, a clustering strategy was also developed and applied to define sub-regions that prioritize one of the objectives and an intermediary region that provides a balance between objectives. This is a valuable tool in the context of process optimization, aiming at assertive decision-making purposes. As this is a novel concept, the only way to compare it was to draw the entire regions of the benchmark functions and compare them with the methodology result. The benchmark results demonstrated that the proposed method could efficiently portray the Pareto regions. Then, the optimization of a Pressure Swing Adsorption unit was performed using the proposed approach to provide a practical application of the methodology developed here. It was possible to build the Pareto region and its respective sub-regions, where each process performance parameter is prioritized. The results demonstrated that this methodology could be helpful in processes optimization and operation. It provides more flexibility and more profound knowledge of the system under evaluation.

1. Introduction

In the optimization of chemical processes, several problems are confronted. These can be formulated from a single-objective and multiobjective point of view. The choice of one of these approaches considers the relationship between the problem and the most appropriate solution. In other words, the type of problem will define the most appropriate strategy to solve it. However, a multiobjective approach is the most suitable when dealing with conflicting objectives, and a trade-off needs to be evaluated.

In fact, the multiobjective optimization (MOO) technique allows the determination of possible conditions portrayed in a set of solutions known as Pareto-optimal solutions. The Pareto front contains the evolution of a given goal to the detriment of another one. It is a valuable tool vastly employed in the optimization literature [1,2,3,4].

It has been noticed recently in the single-objective optimization literature a tendency towards developing algorithms that provide the optimal point and an optimal region. For instance, [5] presented a methodology for the feasible operation region map of the optimal point for single-objective unconstrained problems. These optimal operating regions are named by the authors as feasible operation region (FOR). The method presented in the referred work was extended by [6] to the context of constrained single-objective optimization. The FORs are demonstrated to be a useful tool for system optimization, which makes use of already available results to transform them into decision-making tools. However, in the context of multiobjective optimization problems, there is still a lack of methodology for building the feasible operation region from Pareto fronts for both constrained and unconstrained cases. Therefore, it is an open issue in a multiobjective scenario.

The construction of the feasible operation region can be made possible by the use of metaheuristic optimization methods, such as Evolutionary Algorithms (EAs), Genetic Algorithms (GA), and Particle Swarm Optimization (PSO). As a result, these optimization methods offer an extensive database that contains essential information about the system behavior. This information is readily available to be used to trace a feasible operation region [7]. In addition, it is a waste of knowledge and computational effort not to take advantage of the potential that the data population can offer for the operation of the process.

There is increasing complexity in the problems dealt with in several engineering fields. Most of this complexity comes from Industry 4.0 revolution. Following the increase of complexity, systems optimization has difficulty keeping track of several different demands. This requires new optimization strategies that can cope with these changes. This issue has been receiving attention in recent publications. For example, [8,9] investigated a complex multiobjective problem of 3D redeploying in indoor connected nodes in the Internet of Things collection networks. The authors point out the necessity of developing modified optimization algorithms to solve the above-mentioned complex problems in the referred works. The authors of [7] study the optimization of a complex multifunctional chemical problem. In the referred work, the authors also pointed the necessity of developing a new optimization algorithm that can deal with the peculiarities and complexities of such a system.

On the other hand, there is an increasing concern in the literature to improve Pareto optimization strategies. For instance, in [10], a method was developed to determine Pareto fronts’ uncertainty by defining the optimal upper and lower limits in the objective space. The Pareto uncertainty evaluation proposed in the referred work was applied to address an optimization problem for land use allocation (multiobjective problem with constraints). Their approach uses optimal solutions combined with the outdoor sampling method to propagate an uncertainty of the spatial input data. In [11], a method is presented for evaluating the Pareto front’s uncertainty based on conditional simulations of Gaussian processes. The referred work considers the Kriging metamodels as a tool to estimate the entire Pareto front and quantify the uncertainties at any stage of Kriging-based multiobjective optimization algorithms. Another further effort in calculating the Pareto front uncertainty is described by [12]. Even though uncertainty evaluation is an essential step in this field, it is limited to the perspective of the prior sources of uncertainties.

To date, there is a lack of works that deal with mapping all possible sets of decision variable points that will lead to an optimal value equivalent to the Pareto front’s non-dominant set. This posterior evaluation of the optimization results is proposed in the present work. It is important to highlight that when one performs a traditional Pareto optimization, the algorithm is only concerned with identifying the Pareto front. If a population-based method is applied, the population resulting from this procedure is discarded. However, there is no significant extra computation effort in performing simple mathematical operations to analyze the readily available results. Thus, this work proposes a strategy that analyzes these results to transform them into a feasible operating region (FOR) instead of a single Pareto front. Here lies the main contribution of this work, proposing a methodology that, through a simple procedure, makes it possible to analyze the optimization results of a meta-heuristic algorithm and transform them into useful tools. This concept of FORs is based on the search within the optimization results to identify similar operating conditions that will lead to an equivalent optimal Pareto’s point with a statistical degree of confidence. This provides more flexibility to decision-making and systems operation.

This flexibility is helpful in terms of the practical application of the optimization results. Instead of operating a system using a single line and its corresponding conditions, one can determine some other possibilities that the system can operate. Furthermore, through the proposed method, it is possible to identify the areas and set of conditions that will prioritize a given objective.

In this scenario, the present work aims to develop a concise and systematic methodology to assess feasible operation regions for multiobjective optimization problems. A Fisher–Snedecor test is then deduced to identify FOR from the constrained multiobjective optimization results obtained through meta-heuristic algorithms. The deductions are based on an algebraic point of view to avoid any statistical assumption that would oversimplify the issue. Furthermore, a clustering technique is introduced in this context to leverage the possibilities of the new Pareto regions. Thus, the region is clustered into subregions that offer a priority to a given objective. This can be a helpful tool in, for example, processes operation. A series of four benchmark tests are conducted to evaluate the efficiency of the proposal. A final case study is presented where the method is applied to optimize a complex chemical process.

Given the present scenario, the main contribution of this work consists of developing a concise and systematic methodology to estimate the shape of the Pareto-optimal region, drawing the so-called FORs in a Pareto context. The method proposed here enriches the discussion in the field, delivering a more comprehensive multiobjective optimization approach.

2. Materials and Methods

2.1. Optimization Algorithm

It is sought in the present work to develop a methodology that uses an already available optimization result without the necessity of additional evaluations of the objective function. Therefore, the population-based metaheuristic methods are the most suitable approach to be used in this work. Among these kinds of algorithms, the PSOs are highlighted as being an efficient option.

In terms of classification, PSO is considered an easy-to-implement optimization method, applicable to general problems, and presents better results when compared to GA in several applications [5]. From the point of view of convergence, the PSO has a linear relationship between population size and solution time, while GA requires much greater computational effort to solve similar problems [5,13].

The work [6] presents a novel Constrained Sliding Particle Swarm Optimizer (CSPSO). In the referred work, the CSPSO behavior was evaluated in large-size problems and compared with alternative methods. The authors have demonstrated the efficiency of this algorithm to solve constrained single-objective problems while providing their feasible operating regions (FOR)s. Therefore, this work uses the CSPSO as an optimization strategy to address the evaluated problem, focusing on development of a strategy to build FORs in multiobjective optimizations procedures. The CSPSO was used here due to its advantage of dealing with multimodal problems. For more details regarding this algorithm, please refer to [6].

2.2. Mapping Pareto Region

The concept of feasible operation regions analysis of an optimization result allows deepened insights about the optimal point, which instead of a fixed point becomes an area counterpart. In the context of process systems engineering, the FORs provide flexibility to the process operation. Furthermore, as will be presented, they allow data recycling, using the large amount of data generated by an optimization procedure.

To achieve a reliable FORs, a statistical evaluation of the population given by the optimization is necessary. In many works, the Fisher–Snedecor criterion is used for similar purposes [5,6,14,15,16]. However, it is not available in the literature a Fisher–Snedecor test that considers constrained multiobjective optimization. Therefore, the present work proposes to deduce this test algebraically. As we are dealing with constrained multiobjective optimization problems, Lagrangian functions were considered as described in Equation (1):

(1)$L_i\left(\mathit{\theta },\mathit{\lambda }\right)=fo{b}_i\left(\mathit{\theta }\right)+{\mathit{\lambda }}^{\top }\mathit{f}\left(\mathit{\theta }\right)$

where $\mathit{\lambda }$ is known as a vector of Lagrange multipliers; $fo{b}_i$ is the i-th objective function, $\mathit{\theta }$ is a vector of decision variables and $\mathit{f}$ is the vector of constraints.

Considering a theoretical Pareto curve where a multiobjective optimization will lead to the Pareto global front, in the extremes of that front, it is possible to find the global optima for each objective: $L_A\left({\mathit{\theta }}_{ }^g,{\mathit{\lambda }}^g\right)$and $L_B\left({\mathit{\theta }}_{ }^g,{\mathit{\lambda }}^g\right)$ with their correspondent set of decision variables ${\mathit{\theta }}_A^g$ and ${\mathit{\theta }}_B^g$. Please note that for the sake of the visualization, the deductions are made from using two objectives, but these deductions are valid for n-objectives.

In Figure 1, ${\in }_A$ and ${\in }_B$ correspond to the difference between the global optima for a given goal and the optimal points found. Now consider a particular set of decision variables, ${\mathit{\theta }}_k$, found during optimization by the CSPSO that refers to the set of dominated points, which will have axial coordinates $L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)$ = [$L_A\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)$ $L_B\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right) ]$ $\forall k=1, 2, 3, \dots$, n_k where n_k is the number of population points.

From the previous concepts, it is possible to compute the scaled distance between a point located in the population and an optimal point in the Pareto front, which is here represented by $e_i^{*}$. Therefore, the scaled squared error between those optimal points is given by:

(2)$e_{ij}^{*}{}^2=\frac{{\in }_{ij}^{*}{}^2}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)}={\displaystyle \sum }_{k=1}^{n_k}\frac{{\left(L_i({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})-L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)\right)}^2}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)},\forall i=A,B,$

where $V_A$and $V_B$are the respective variances of $L_A({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})$ and $L_B({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})$ values. Then, the error can be generalized along the optimized Pareto front as:

(3)$e_i^2={\displaystyle \sum }_{j=1}^{n_j}{e}_{ij}^{*}{}^2=\frac{1}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)}{\displaystyle \sum }_{j=1}^{n_j}{\displaystyle \sum }_{k=1}^{n_k}{\left(L_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)-L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)\right)}^2$

where $i$ is the corresponding objective, $A$ or $B$ and $n_j$ is the overall number of optimal points in the Pareto front. Considering now that the variance of the optimal points found can be calculated as a function of the population deviation of dominated points as:

(4)$V_{ij}\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)=\frac{1}{n_k}\left({\displaystyle \sum }_{k=1}^{n_k}{\left(L_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)-L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)\right)}^2\right)$

While the variance $V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)$, computed about the theoretical Pareto front and the obtained Pareto front, is given by the expression that follows:

(5)$V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)=\frac{1}{n_j}\left({\displaystyle \sum }_{k=1}^{n_j}{\left(L_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)-L_i\left({\mathit{\theta }}^g,{\mathit{\lambda }}^g\right)\right)}^2\right)$

Hence the deviation can be presented as:

(6)$e_i^2={\displaystyle \sum }_{j=1}^{n_j}{e}_{ij}^{*}{}^2={\displaystyle \sum }_{j=1}^{n_j}\frac{\left(n_k\right)V_{ij}\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)}$

Assuming that $e$ follows a Gaussian distribution with null mean, then the $e_i^2$ will follow a chi-squared distribution, ${\chi }^2$, having $n_k-\left(n_j+n_{\theta }-1\right)$ degrees of freedom, with $n_{\theta }$ equal to the number of decision variables, thus:

(7)$e_i^2 \to {\chi }^2\left(n_k-n_j-n_{\theta }+1\right)$

Please note that as it is an optimization problem, the decision variables will tend to be randomly distributed around the optimal point. Thus, the null mean assumption is reasonable in this context. Now, considering the difference between the theoretical optimal points and a subsequent point from the cloud found during the optimization as:

(8)$b_{ij}^2={\displaystyle \sum }_{k=1}^{n_k}\frac{{\left(L_i\left({\mathit{\theta }}^g,{\mathit{\lambda }}^g\right)-L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)\right)}^2}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)}$

Extending this towards all points obtained in each objective function evaluation made by the optimization, the scaled difference can be represented as:

(9)$b_i^2={\displaystyle \sum }_{k=1}^{n_k}{h}_{ij}^{*}{}^2={\displaystyle \sum }_{k=1}^{n_k}\frac{\left(n_k\right){\widehat{V}}_{ij}\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)}{V_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)}$

and g is calculated as:

(10)${\widehat{V}}_{ij}\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)=\frac{1}{n_k}\left({\displaystyle \sum }_{k=1}^{n_k}{\left(L_i\left({\mathit{\theta }}^g,{\mathit{\lambda }}^g\right)-L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)\right)}^2\right)$

where$, n_k$ is the number of total evaluations of the objective function made in the optimization, $n_g$ is the overall number of optimal points in the theoretical Pareto front. Therefore, $h$ follows a chi-squared distribution as well, having $n_k-\left(n_{\theta }+n_g-1\right)$ degrees of freedom, and is approximated as:

(11)$b_i^2 \to {\chi }^2\left(n_k-n_{\theta }-n_g+1\right)$

where $n_{\theta }$ is the number of optimized operating variables. Thus:

(12)$b_i^2-e_i^2\to {\chi }^2\left(n_j-n_g\right)$

Geometrically, if $b_i^2$ is a chi-squared distribution, then ${\widehat{b}}_i^2$ is also a chi-squered distribution proportional to $b_i^2$. Moreover, it is known that there is a geometrical relationship between $b$ and $h$. Thus, the approximation of the distances as below leads to a ratio of two chi-squared distributions:

(13)$\frac{h_i}{{\widehat{b}}_i}=\frac{b_i}{h_i}=\frac{b_i^2}{e_i^2-{\widehat{b}}_i{}^2}$

From the Pythagorean theorem $e_i^2-{\widehat{b}}_i{}^2$is equal to $h_i^2.$ Furthermore, $h_i$ and $b_i$ are orthogonal to each other, as seen in Figure 1. Since the vectors are independent, they will also be independent distributions, so the proportion is equivalent to a Fisher–Snedecor distribution, as:

(14)$\frac{\frac{b_i^2}{n_k-n_{\theta }-n_g+1}}{\frac{e_i^2-{\widehat{b}}_i^2}{n_j-n_g}}\to F_{\alpha }\left(n_k-n_{\theta }-n_g+1,n_j-n_g\right)$

Finally, Equation (6) can be expressed for the Equation (14) as:

(15)${\displaystyle \sum _{k=1}^{n_k}\frac{(n_k){\widehat{V}}_{ij}({\mathit{\theta }}_k,{\mathit{\lambda }}_k)}{V_i({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})}}\to {\displaystyle \sum _{j=1}^{n_j}\frac{(n_k)V_{ij}({\mathit{\theta }}_k,{\mathit{\lambda }}_k)}{V_i({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})}}-{\displaystyle \sum _{k=1}^{n_k}\frac{(n_k){\widehat{V}}_{ij}({\mathit{\theta }}_k,{\mathit{\lambda }}_k)}{V_i({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*})}}\frac{n_k-n_{\theta }-n_g+1}{n_j-n_g}{F}_{\alpha }(n_k-n_{\theta }-n_g+1,n_j-n_g)$

where $\alpha$ is the confidence level of the Fisher–Snedecor test. Deriving Equation (15) with respect to $L_i\left({\mathit{\theta }}_k,{\mathit{\lambda }}_k\right)$ will provide the level of approximation of a given point computed by the optimizer to the optimal point. Furthermore, assuming that the optimal value is calculated by the optimizer, $L_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)$, will tend to the theoretical value $L_i\left({\mathit{\theta }}^g,{\mathit{\lambda }}^g\right)$, and this tendency can be represented by a factor $\beta$, according to Equation (16):

(16)$L_i({\mathit{\theta }}_k,{\mathit{\lambda }}_k)\le L_i\left({\mathit{\theta }}_j^{*},{\mathit{\lambda }}_j^{*}\right)+\frac{\beta (n_k-n_{\theta }-n_g+1)}{n_j-n_g}{F}_{\alpha }\left(n_k-n_{\theta }-n_g+1,n_j-n_g\right)$

Hence, the extended Fisher–Snedecor test is reached to a multiobjective optimization problem applied to the Pareto front. This test can be used to evaluate a population obtained in a multiobjective optimization, which will construct the FOR of the optimal Pareto front. Thus, the Pareto set can be represented as a curve and region, a more flexible concept.

2.3. Clustering

Once the optimal Pareto region is built, it can be divided into subregions, using a clustering method to map the Pareto region. For example, this tool is interesting for real processes. It provides relevant information about optimal process conditions with greater flexibility for system operating parameters.

Figure 2 shows a schematic diagram of the approach proposed here. Figure 2a shows the optimal Pareto region defined according to the Fisher–Snedecor criterion. In turn, in Figure 2b, the optimal region is divided into sub-regions with different commitments regarding objective functions, using a clustering technique that identifies three centroids within the FORs. Thus, as they are conflicting goals, the Pareto region will present sub-regions, varying priorities for the objectives. Finally, Figure 2c illustrates how these sub-regions will affect the FORs of the set of decision variables.

There are two objectives in the presented optimization problem for this optimization problem, and the objective function is a continuous function. From the decision-maker perspective, we would have three options, an action that prioritizes one objective to the detriment of the other, the inverse action, or an action that would bring a balance between the two objectives. Therefore, it was used three clusters represent this chain of possible actions.

Therefore, to define the clusters of the optimal Pareto region, a metric must be selected. Various metrics may be suitable, such as the unweighted pair group method using centroids. This method consists of searching subgroups within a region through the definition of a center of mass. Then the Euclidean distance is used to define which points have the smallest distance from the center. The point closer to a center of mass will compose a sub-region. Figure 3 presents a simplified flowchart that exemplifies a metric that can define the clusters through the centroid.

3. Results

A series of benchmark tests were performed to assess the consistency of the feasible operational regions proposed in this work. The first step was to define the CSPSO optimization parameters based on [6], presented in Table 1.

The confidence level is associated with the desired region dimension that the user wants for results. The higher the level, the shorter will be the FOR.

For each benchmark, this algorithm was executed to solve the related optimization problem. The final population, provided by the CSPSO, was ranked by the test here proposed, and the FORs were drawn. Furthermore, as there is no alternative approach to perform the proposal of this work, Monte Carlo regions were drawn for the benchmarks and compared with the proposed method. Please note that the Monte Carlo region is computationally hard to obtain. Therefore, it was not performed for the case study, as it is by itself a computationally heavy problem.

3.1. Benchmark 1—SRN

The first benchmark function tested is the SRN, a convex, nonlinear, constrained, two-parameter function [17,18]. The SRN benchmark consists of two objective functions described as follows:

(17)$Minimize F=\left(f_1\left(x\right), f_2\left(x\right)\right),$

(18)$\begin{array}{c}{f}_1\left(x_1,x_2\right)={\left(x_1-2\right)}^2+{\left(x_2-1\right)}^2+2\end{array}$

(19)$f_2\left(x_1,x_2\right)=9x_1-{\left(x_2-2\right)}^2$

Subject to

(20)$g_1\left(x_1,x_2\right)=x_1^2+x_2^2\le 225$

(21)$g_2\left(x_1,x_2\right)=x_1-3x_2+10\le 0$

(22)$-20\le x_1,x_2\le 20$

This benchmark problem is composed of two functions, one with a monotonously decreasing behavior with the decrease of $x_1$ (and the increase of the absolute value of $x_2$) and another with a smooth unimodal behavior. The interaction between these two objective functions favors the construction of a predominantly linear Pareto front, with some nonlinearities at its extremes.

Figure 4 shows the Pareto front obtained for the optimization problem; the optimal Pareto region is divided into sub-regions by clustering. Furthermore, in Figure 4, the whole feasible region of this objective function is depicted through Monte Carlo. As it is possible to see, the Pareto region drawn by the proposed methodology covers and enters, the left limit of the MC area. This overlapping evidences that the proposed method can efficiently identify the Pareto limit. Furthermore, the identified Pareto region is within the MC region, indicating that the proposed methodology can approximate the Pareto region.

As a consequence of the proposed methodology, it is possible to draw the feasible operation regions of the decision variables, portrayed in Figure 5, together with its respective sub-regions. This figure shows the relationships between the decision variables and their influence in the sub-regions of the Pareto region. It is interesting to note the function nonlinearities expressed in terms of FOR. Furthermore, based on these FOR, it is possible to identify the set of decision variables that will lead to an equivalent value of the objective function or an equivalent value of sub-region.

3.2. Benchmark 2—TNK

The TNK function was proposed by [19,20,21]. This test function presents a discontinuous Pareto front as a solution to the optimization problem, which is nonlinear and constrained. Furthermore, it does not present a linearly distributed solution along the Pareto front, which is an attractive test for evaluating the proposed methodology. The objective functions are described as follows:

(23)$Minimize F=\left(f_1\left(x_1\right), f_2\left(x_i\right)\right)$

(24)$\begin{array}{c}{f}_1\left(x_1\right)=x_1 \end{array}$

(25)$f_2\left(x_2\right)=x_2$

where

(26)$g_1\left(x_1,x_2\right)=x_1^2+x_2^2-1-0.1\cos \left(16\arctan \left(\frac{x_1}{x_2}\right)\right)\ge 0$

(27)$g_2\left(x_1,x_2\right)={\left(x_1-0.5\right)}^2+{\left(x_2-0.5\right)}^2\le 0.5$

(28)$0\le x_1,x_2\le \pi$

Figure 6 depicts the Pareto region obtained after solving this optimization problem through the proposed methodology. As in the previous case, the optimal Pareto region is divided into sub-regions by clustering. As a test to evaluate the results obtained, the feasible area of this objective function was also portrayed by Monte Carlo simulations. As the results obtained in the previous case, the Pareto region drawn by the proposed methodology covers the entire left limit of the MC area, overlapping part of the left region, always within the MC region. The observed overlapping corroborates the previously obtained evidence; the proposed method efficiently approximates the Pareto region. As this benchmark is an irregular, constrained function, these observations are important in validating the methodology proposed here.

These results are again transported to the decision variables context. Figure 7 depicts the variables’ FORs together with its respective sub-regions. Similarly, the function nonlinearities are evident in the graphic, which is a crucial conclusion. The method is capable of reproducing the nonlinearities into the feasible operation regions. As aforementioned, through the FOR, it is possible to identify the set of decision variables that will lead to an equivalent value of the objective function or an equivalent value of sub-region.

3.3. Benchmark 3—ZDT3

For this work, peculiarities of the ZDT family (Zitzler–Deb–Thiele, ZDT3) were also tested [22,23,24]. This function is widely used to evaluate the performance of multiobjective optimization methods. Selecting the ZDT3 is interesting because it adds a discreteness feature to the front, which is unusual.

The ZDT3 Pareto-optimal front consists of several noncontiguous and nonlinear convex parts. The introduction of a sine function in this objective function causes discontinuities in the Pareto-optimal front. The objective functions are described as follows:

(29)$Minimize F=\left(f_1\left(x_1\right), f_2\left(x_1,x\right)\right)$

(30)$\begin{array}{c}{f}_1\left(x_1\right)=x_1 \end{array}$

(31)$f_2\left(x_1,x\right)=g\left(x\right)h\left(f_1\left(x_1\right),g\left(x\right)\right)$

(32)$g\left(x\right)=1+\frac{9}{\left(n-1\right)}{\displaystyle \sum }_{i=2}^n{x}_i$

(33)$h\left(f_1\left(x_1\right),g\left(x\right)\right)=1-\sqrt{\frac{f_1\left(x_1\right)}{g\left(x\right)} }-\left(\frac{f_1\left(x_1\right)}{g\left(x\right)}\right)\sin \left(10\pi f_1\left(x_1\right)\right)$

where

(34)$x=\left[x_2,\dots ,n\right]$

(35)$0\le x_i\le 1$

(36)$1 \le n \le 30$

Figure 8 presents this peculiar Pareto region that is obtained after solving the ZDT3 optimization problem, with $n=2$. Following the previous evaluation, the complete area of this objective function was identified using Monte Carlo simulations. As the results obtained in the last cases, the Pareto region drawn by the proposed methodology is overlapping the Monte Carlo area. Even in this complex subject, the proposed method can provide an approximation of the Pareto regions. This observation is a piece of more robust evidence that the methodology can identify the Pareto region. Moreover, the behavior of this function allows the visualization in Figure 8 of how the Pareto is located along the objective function feasible region. It is also possible to see how the proposed methodology can track the Pareto approximating a corresponding surface of FORs.

Finally, the Pareto region is used to build the decision variables feasible operation region. Figure 9 presents the variables’ FORs together with its respective sub-regions. Interestingly, this problem results in a series of discrete FOR. The peculiar behavior of this function is highlighted in the graphic. The method is capable of reproducing the nonlinearities into the feasible operation region. As aforementioned, through the FOR, it is possible to identify the set of decision variables that will lead to an equivalent value of the objective function or an equivalent value of sub-region.

3.4. Benchmark 4—Modified OffsetValleys

The OffsetValleys function is based on a test function proposed in [25]. This benchmark is based on the expansion of terms of a Gaussian basis function. In this work, the formulation was extended to represent a multiobjective problem with two decision variables, namely:

(37)$$\begin{gathered} Minimize F=\left(f_1\left(x_1,x_2\right), f_2\left(x_1,x_2\right)\right) \\ \mathrm{With}, \end{gathered}$$

(38)$$\begin{gathered} f_1\left(x_1,x_2\right)=-0.6e^{-\frac{{\left(x_1-1\right)}^2+{\left(x_2-1\right)}^2}{0.18}}-0.75e^{-\frac{{\left(x_1-1\right)}^2+{\left(x_2-3.5\right)}^2}{0.32}} \\ -0.9e^{-\frac{{\left(x_1-3\right)}^2+{\left(x_2-1\right)}^2}{2}}-1.2e^{-\frac{{\left(x_1-4\right)}^2+{\left(x_2-4\right)}^2}{0.32}} \end{gathered}$$

(39)$$\begin{gathered} f_2\left(x_1,x_2\right)=-0.6e^{-\frac{{\left(x_1-1.2\right)}^2+{\left(x_2-1.2\right)}^2}{0.18}}-0.75e^{-\frac{{\left(x_1-1.2\right)}^2+{\left(x_2-3.7\right)}^2}{0.32}} \\ -0.9e^{-\frac{{\left(x_1-3.5\right)}^2+{\left(x_2-1.3\right)}^2}{2}}-1.2e^{-\frac{{\left(x_1-3.8\right)}^2+{\left(x_2-3.8\right)}^2}{0.32}} \end{gathered}$$

where

(40)$0\le x_1,x_2\le 5$

Each of the four exponential clusters produces a domain of attraction accountable for generating a series of local and global Pareto-optimal fronts, as seen in Figure 10. This problem places a new challenge for the methodology proposed here, which is coping with a multimodal and multiobjective constrained problem. As is possible to see in Figure 10, the strategy developed here was able to identify each local Pareto and subdivide these Pareto into the sub-regions.

To evaluate the efficiency of the proposed methodology solving this problem Figure 11 is presented. In Figure 11a it is shown the surface map of the analytical solution of each function involved in this problem. Therefore, it is possible to analyze where the feasible operation region are located concerning the surface of the objective function. This graphic shows that the feasible operation regions drawn by the deduced test are within the minimum of each objective function. Moreover, it is possible to see that the sub-regions are indeed prioritizing one objective.

In Figure 11b, the 2D and 3D maps of Figure 11 provide a complete picture of the Modified OffsetValley function topography. The 2D map is a projection of the 3D topography in a perpendicular plane. It is essential to remind that the Pareto front is defined by one trade-off between two or more objectives. In this way, Figure 11a provides an interesting analysis through the contour’s curves. The “puddles”, minimum local areas, of the functions work like a “decoy” to the objective function, and when the optimizer tries to go to the center of a “puddle” of one objective, it will dissipate the other. This effect is apparent in Figure 11. For instance, the sub-region 1 is always closer to the minimum of the objective function 1 being out of the minimum of the objective function 2. In contrast, it is possible to see the opposite effect in region 3. Thus, region 2 is located in the middle way between the objective 1 local minimum and the local minimum of objective function 2. Hence, the three sub-regions together form a complete picture of the “puddles” of both objective functions. From these analyses, it is possible to confirm that the methodology proposed here can provide a reasonable estimation of the Pareto extension.

3.5. Performance Comparison

A metrics to evaluate the performance of an algorithm is essential to compare it. Several efficiency metrics for multiobjective algorithms have been proposed in the literature [22,26,27]. In this work, five evaluation metrics are used to evaluate the CSPSO performance, namely: generational distance (GD) [28], generational distance plus (GD+) [29], inverted generational distance (IGD) [30], inverted generational distance plus (IGD+) [29], and hypervolume (HV) [31]. A traditional Genetic Algorithm was also applied to solve the benchmark problems, serving as a reference to compare the CSPSO performance. The GA was applied with 100 individuals in the population over 100 generations. While the CSPSO with 100 particles and 100 iterations. The metrics GD, GD+, IGD, IGD+ depend on the known Pareto front value and the solutions obtained for each algorithm. Table 2 provides the performance metrics for both CSPSO and GA, computed for each benchmark function tested. That table also presents the number of points on the Pareto front obtained by the different optimization algorithms. This last metric is an essential reference for the decision-making step. As it provides an overview of the number of possible conditions the system can be operated with. The lower the feasible operating region population, the lower is the flexibility for the system operation.

In general, the CSPSO obtained better IGD and IGD+ metrics. In contrast, the others did not show significant differences, except for the GD and GD+ in the first benchmark. However, the number of points present on the Pareto front shows excellent filling compared to the solution obtained from GA. Furthermore, it is important to highlight that the GA presented a significantly higher computational effort. When compared in terms of computing time, the GA gives a fivefold factor higher than the CSPSO.

3.6. A Case Study—Multiobjective Optimization of a Pressure Swing Adsorption Process

The Fischer–Tropsch process is an industrial process for producing liquid hydrocarbons, such as gasoline, kerosene or diesel, from synthesis gas. As these reactions take place in the presence of a catalyst, they must be protected from impurities that can compromise the reactions yield due to catalytic deactivation. Furthermore, the H₂/CO ratio must be adjusted between 2.0 and 2.4 to ensure that carbon is not deposited on the catalyst, deactivating it. To ensure that Fischer–Tropsch process specifications are met, an additional separation process is necessary to purify its feed stream. Usually, a pressure swing adsorption (PSA) unit can be employed for this goal. Therefore, purifying the synthesis gas by removing impurities and adjusting the unit load.

The pressure swing adsorption process was used as a case study in this work. This process dynamic behavior is described by mass, energy, and momentum balances and constitutive relationships associated with adsorption equilibrium and kinetics, besides unit boundary conditions. The process does not reach a steady state; instead, it has a cyclic steady state. This is a very complex dynamic system, which turns out to be an excellent practical case to evaluate the method proposed here. The work [32] presented the parameter identification for the syngas purification in a PSA, validating the mathematical model against experimental data. Please refer to this work if the reader is interested in more details about the process, phenomenological model, and the system.

This model is used here as a case study, where a multiobjective optimization problem is proposed to exemplify the existing trade-off between the CO₂ Purity ($Pu{r}_{{\mathrm{CO}}_2}$) and the H₂/CO productivity ($P{r}_{{\mathrm{H}}_2+\mathrm{CO}}$). The vector of decision variables $\mathit{\theta }$ is composed of: the duration of three steps ($t_{feed}, t_{rinse}, t_{purge}$), operating pressures ($P_{high} \mathrm{and} P_{low}$), rinse and purge flow rates ($Q_{rinse}, Q_{purge}$), and inlet temperature ($T_{inlet}$). The search region limits provided to the CSPSO optimizer are given in Table 3.

The multiobjective optimization problem was formulated as follows:

(41)$\{\begin{array}{c}\underset{\mathit{\theta }}{\mathrm{Max}}Pu{r}_{{\mathrm{CO}}_2}\left(\mathit{\theta },\mathit{p}\right) \\ \underset{\mathit{\theta }}{\mathrm{Max} } P{r}_{{\mathrm{H}}_2+\mathrm{CO}}\left(\mathit{\theta },\mathit{p}\right) \end{array}$

Subject to:

(42)$\mathit{\theta } \in [{\mathit{\theta }}_{min}, {\mathit{\theta }}_{max}]$

(43)$2.2\le \frac{{\mathrm{H}}_2}{\mathrm{CO}}\left(\mathit{\theta };\mathit{p}\right)\le 2.3$

(44)$Pu{r}_{{\mathrm{CO}}_2}\left(\mathit{\theta };\mathit{p}\right)\ge 84\%$

(45)$Re{c}_{{\mathrm{CO}}_2}\left(\mathit{\theta };\mathit{p}\right)\ge 95\%$

where:

(46)$\begin{array}{l}P{r}_{{\mathrm{H}}_2+\mathrm{CO}} \left[{\mathrm{mol}}_{{\mathrm{H}}_2+\mathrm{CO}}{ \mathrm{kg}}^{-1}{ \mathrm{h}}^{-1}\right]\\ =\frac{{\int }_0^{t_{feed}}{C}_{{\mathrm{H}}_2+\mathrm{CO}}{u}_0{|}_{z=L}dt+{\int }_0^{t_{rinse}}{C}_{{\mathrm{H}}_2+\mathrm{CO}}{u}_0{|}_{z=L}dt-{\int }_0^{t_{purge}}{C}_{{\mathrm{H}}_2+\mathrm{CO}}{u}_0{|}_{z=L}dt}{t_{cycle}\times mass dry adsorbent}\end{array}$

(47)$Pu{r}_{{\mathrm{CO}}_2}=\frac{{\int }_0^{t_{blowdown}}{C}_{{\mathrm{CO}}_2}{u}_0{|}_{z=0}dt+{\int }_0^{tpurge}{C}_{{\mathrm{CO}}_2}{u}_0{|}_{z=0}dt}{{\displaystyle \sum _{i=1}^n}[ {\int }_0^{t_{blowdown}}{C}_i{u}_0{|}_{z=0}dt+{\int }_0^{t_{purge}}{C}_i{u}_0{|}_{z=0}dt ]}$

(48)$\frac{{\mathrm{H}}_2}{\mathrm{CO}}=\frac{{\int }_0^{t_{feed}}{C}_{{\mathrm{H}}_2}{u}_0{|}_{z=L}+{\int }_0^{t_{rinse}}{C}_{{\mathrm{H}}_2}{u}_0{|}_{z=L}}{{\int }_0^{t_{feed}}{C}_{\mathrm{CO}}{u}_0{|}_{z=L}+{\int }_0^{t_{rinse}}{C}_{\mathrm{CO}}{u}_0{|}_{z=L}}$

(49)$\begin{array}{l}Re{c}_{\mathrm{C}{\mathrm{O}}_2}\\ =\frac{{\displaystyle {\int }_0^{t_{blowdown}}{C}_{\mathrm{C}{\mathrm{O}}_2}{u}_0|{}_{z=0}dt+{\displaystyle {\int }_0^{t_{purge}}{C}_{\mathrm{C}{\mathrm{O}}_2}{u}_0|{}_{z=0}dt-{\displaystyle {\int }_0^{t_{rinse}}{C}_{\mathrm{C}{\mathrm{O}}_2}{u}_0|{}_{z=0}dt}}}}{{\displaystyle {\int }_0^{t_{press}}{C}_{\mathrm{C}{\mathrm{O}}_2}{u}_0|{}_{z=0}dt+{\displaystyle {\int }_0^{t_{feed}}{C}_{\mathrm{C}{\mathrm{O}}_2}{u}_0|{}_{z=0}dt}}}\end{array}$

Additionally θ is the vector of decision variables, ${\mathit{\theta }}_{min}$ and ${\mathit{\theta }}_{max}$ are the lower and upper limits of the decision variables, p is the parameter set, $u_0$ is the superficial velocity, $C_i$ is the concentration of chemical species i, z is the axial position, L is the bed length,$t_{press}, t_{feed}, t_{rinse}, t_{blowdown},t_{purge}$ correspond to the time of the pressurization, adsorption, rinse, blowdown, and purge operation steps.

The optimization problem was solved using 100 iterations and 100 particles to build the Pareto region shown in Figure 12. As it is possible to see the result of this optimization problem through the methodology proposed here is an extended Pareto front of this chemical process, where all the points are equivalently leading to an optimal operation. Furthermore, the clustering technique split the Pareto region into three sub-regions, each providing a different commitment to the objectives. This concept of sub-region is a valuable tool in decision-making. For instance, one can quickly know where to operate the process to prioritize productivity, Region 1 in Figure 12; to provide a balance between productivity and purity, Region 2; or where prioritize the purity, Region 3.

Following the methodology proposed here, it is possible to build the feasible operation region of the decision variables. Now in an n-dimensional scenario. The complex relationships between the operating variables of this process and how they relate to the prioritized objectives are then mapped in Figure 13. An interesting analysis that can be extracted from Figure 13 is that all the sub-regions are mixed except the regions concerning the duration of the rinse step. This provides extra insights into the process operation. One can see that the rinse duration is the only variable that can be manipulated to move from a sub-region to another.

Hence, the methodology proposed in this work provides an extra level of understanding and flexibility to the process’s optimization context. It allows visualizing how the operating variables are related to possible ways to operate the process.

4. Conclusions

This work contributes to the optimization literature proposing a new concept of feasible operation regions estimation for multiobjective optimization problems. The methodology is based on the algebraic deductions of a novel Fisher–Snedecor test. Through this test, it is possible to evaluate the population of points generated by metaheuristic optimizations. Hence, it is possible to estimate a feasible operation region that extends the optimal Pareto front to an optimal Pareto region.

Multiobjective optimization problems have two essential steps: (i) identifying the optimal solution, which is the Pareto set, and (ii) the decision making based on the Pareto solution, as it typically portrays a conflict between the objectives of the problem. The methodology presented in this work contributes to the two steps. In the first part, the novel concept of Pareto regions is proposed. In the second step, the methodology contributes to the clustering of the optimal Pareto region, dividing it into sub-regions with different commitments to the objectives in evaluation. Gliding between the different sections of the Pareto optimal region allows one to move between the commitments. The case study demonstrated that the proposed method is an important tool for decision-making.

A series of four benchmark tests with distinct characteristics were used, each with peculiarities that provide a consistent and robust evaluation of the proposed methodology. For all benchmarks tested, the proposed methodology provided a reasonable estimation of the Pareto region. These results were compared against the whole objective function feasible regions drawn by Monte Carlo.

Finally, the optimization of a Pressure Swing Adsorption unit was used as a case study. Applying the proposed methodology to this problem allowed us to verify the methodology in optimizing a complex chemical process. This concept was extended to build the decision variables’ feasible operating regions, providing flexibility to the decision-making based on optimization results. The integration of these concepts into the decision-making step should be further investigated in future works. Furthermore, this work may open new possibilities in the field of real-time optimization, such as employing the idea of multiple FORs in the real-time context. As it is a novel contribution in multiobjective optimization, it should have several potential implications in future developments in this field.

Hence, the tests performed demonstrated that the method presented in this work could map the Pareto optimal front and estimate its feasible operation region. This concept was extended to build the decision variables’ FOR, which can provide flexibility to the decision-making based on optimization results.

Author Contributions

Conceptualization, I.B.R.N., M.A.F.M., D.D.S., C.M.R.; methodology, I.B.R.N., C.M.R.; writing—original draft preparation, I.B.R.N., C.M.R.; writing—review and editing, I.B.R.N., M.A.F.M., C.M.R. and A.M.R.; supervision, I.B.R.N., J.M.L., A.E.R., A.M.R., M.A.F.M. All authors have read and agreed to the published version of the manuscript.

Funding

This work was financially supported by: Project—NORTE-01-0145-FEDER-029384 funded by FEDER funds through NORTE 2020—Programa Operacional Regional do NORTE—and by national funds (PIDDAC) through FCT/MCTES. This work was also financially supported by: Base Funding—UIDB/50020/2020 of the Associate Laboratory LSRE-LCM—funded by national funds through FCT/MCTES (PIDDAC), Capes for its financial support, financial code 001 and FCT—Fundação para a Ciência e Tecnologia under CEEC Institucional program.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Schematic representation of the proposed criteria, where [Forumla omitted. See PDF.] is the triangle height, and [Forumla omitted. See PDF.] is the triangle base, [Forumla omitted. See PDF.] is the base of the iner triangle, [Forumla omitted. See PDF.] is the hypotenuse.

Figure 2. Schematic diagram of the clustering process of the optimal Pareto region. (a) Pareto front; (b) Centroids identification; (c) Clusters.

Figure 3. Simplified flowchart for determining clusters using centroid metric.

Figure 4. Front and Pareto optimal divided into regions by clustering to SRN.

Figure 5. Regions of optimal parameters divided by clustering to SRN.

Figure 6. Front and Pareto optimal region divided into regions by clustering to TNK.

Figure 7. Regions of optimal parameters divided by clustering to TNK.

Figure 8. Pareto optimal front and region divided into regions by clustering to ZDT3.

Figure 9. Regions of optimal parameters divided by clustering to ZDT3.

Figure 10. Four Pareto fronts of the modified OffsetValleys function.

Figure 11. Modified OffsetValleys function 2D topography drawn by the proposed methodology, (a) 2d contour map, (b) 3d landscape.

Figure 12. Pareto region and sub-operating regions.

Figure 13. Feasible operation region for the case study decision variables.

Table 1

CSPSO Parameters.

Number of Particles	200
Number of Iterations	200
$${\mathit{C}}_{\mathbf{1}\mathit{i}}$$	0.5
$${\mathit{C}}_{\mathbf{1}\mathit{f}}$$	2.5
$${\mathit{C}}_{\mathbf{2}\mathit{i}}$$	2.5
$${\mathit{C}}_{\mathbf{2}\mathit{f}}$$	0.5
α (Confidence Level)	0.99

Table 2

GA and CSPSO algorithm performance metrics for the SRN, TNK, ZDT3 benchmarks.

	GD	GD+	IGD	IGD+	HV	Pareto Front Points
SRN
GA	3.0	2.576	0.959	0.756	25,818.831	100
CSPSO	6.843	5.475	0.152	0.133	25,388.842	2196
TNK
GA	0.005	0.003	0.004	0.003	0.291	100
CSPSO	0.005	0.002	0.002	0.001	0.306	602
ZDT3
GA	0.013	0.011	0.015	0.014	0.812	100
CSPSO	0.006	0.002	0.002	0.001	0.781	979

Table 3

The bounds for all decision variables used in multiobjective optimization by Constrained Sliding Particle Swarm Optimizer (CSPSO).

	$${\mathit{t}}_{\mathit{f}\mathit{e}\mathit{e}\mathit{d}}/$$ $$\left(\mathit{s}\right)$$	$${\mathit{t}}_{\mathit{p}\mathit{u}\mathit{r}\mathit{g}\mathit{e}}/$$ $$\left(\mathit{s}\right)$$	$${\mathit{t}}_{\mathit{r}\mathit{i}\mathit{n}\mathit{s}\mathit{e}}/$$ $$\left(\mathit{s}\right)$$	$${\mathit{P}}_{\mathit{h}\mathit{i}\mathit{g}\mathit{h}}/$$ $$\left(\mathit{b}\mathit{a}\mathit{r}\right)$$	$${\mathit{P}}_{\mathit{l}\mathit{o}\mathit{w}}/$$ $$\left(\mathit{b}\mathit{a}\mathit{r}\right)$$	$${\mathit{Q}}_{\mathit{r}\mathit{i}\mathit{n}\mathit{s}\mathit{e}}/$$ $$\left(\mathit{S}\mathit{L}\mathit{P}\mathit{M}\right)$$	$${\mathit{Q}}_{\mathit{p}\mathit{u}\mathit{r}\mathit{g}\mathit{e}}/$$ $$\left(\mathit{S}\mathit{L}\mathit{P}\mathit{M}\right)$$	$${\mathit{T}}_{\mathit{i}\mathit{n}\mathit{l}\mathit{e}\mathit{t}}/$$ $$\left(\mathit{K}\right)$$
Minimum	380	80	187	3.4	0.55	0.425	0.225	304
Maximum	680	110	253	5.0	1.10	0.575	0.345	350

References

1. Yu, C.; Lu, Y.; Chu, J. Multi-objective Optimization with Combination of Particle Swarm and Extremal Optimization for Constrained Engineering Design. WSEAS Trans. Syst. Control; 2012; 7, pp. 129-138.

2. Capra, F.; Gazzani, M.; Joss, L.; Mazzotti, M.; Martelli, E. MO-MCS, a derivative-free algorithm for the multiobjective optimization of adsorption processes. Ind. Eng. Chem. Res.; 2018; 57, pp. 9977-9993. [DOI: https://dx.doi.org/10.1021/acs.iecr.8b00207]

3. Rangaiah, G.P.; Feng, Z.; Hoadley, A.F. Multi-objective optimization applications in chemical process engineering: Tutorial and review. Processes; 2020; 8, 508. [DOI: https://dx.doi.org/10.3390/pr8050508]

4. Alkebsi, K.; Du, W. Surrogate-assisted multi-objective particle swarm optimization for the operation of CO₂ capture using VPSA. Energy; 2021; 224, 120078. [DOI: https://dx.doi.org/10.1016/j.energy.2021.120078]

5. Nogueira, I.B.R.; Martins, M.A.F.; Requião, R.; Oliveira, A.R.; Viena, V.; Koivisto, H.; Rodrigues, A.E.; Loureiro, J.M.; Ribeiro, A.M. Optimization of a True Moving Bed unit and determination of its feasible operating region using a novel Sliding Particle Swarm Optimization. Comput. Ind. Eng.; 2019; 135, pp. 368-381. [DOI: https://dx.doi.org/10.1016/j.cie.2019.06.020]

6. Rebello, C.M.; Martins, M.A.F.; Loureiro, J.M.; Rodrigues, A.E.; Ribeiro, A.M.; Nogueira, I.B.R. From an Optimal Point to an Optimal Region: A Novel Methodology for Optimization of Multimodal Constrained Problems and a Novel Constrained Sliding Particle Swarm Optimization Strategy. Mathematics; 2021; 9, 1808. [DOI: https://dx.doi.org/10.3390/math9151808]

7. Nogueira, I.B.R.; Viena, V.; Rodrigues, A.E.; Loureiro, J.M.; Ribeiro, A.M. Dynamics of a True Moving Bed Reactor: Synthesis of n-Propyl Propionate and an alternative optimization method. Chem. Eng. Process.-Process Intensif.; 2020; 148, 107821. [DOI: https://dx.doi.org/10.1016/j.cep.2020.107821]

8. Mnasri, S.; Nasri, N.; Van Den Bossche, A.; Val, T. 3D indoor redeployment in IoT collection networks: A real prototyping using a hybrid PI-NSGA-III-VF. Proceedings of the 2018 14th International Wireless Communications & Mobile Computing Conference (IWCMC); Limassol, Cyprus, 25–29 June 2018; pp. 780-785. [DOI: https://dx.doi.org/10.1109/IWCMC.2018.8450372]

9. Mnasri, S.; Nasri, N.; Alrashidi, M.; van den Bossche, A.; Val, T. IoT networks 3D deployment using hybrid many-objective optimization algorithms. J. Heuristics; 2020; 26, pp. 663-709. [DOI: https://dx.doi.org/10.1007/s10732-020-09445-x]

10. Hildemann, M.; Verstegen, J.A. Quantifying uncertainty in Pareto fronts arising from spatial data. Environ. Model. Softw.; 2021; 141, 105069. [DOI: https://dx.doi.org/10.1016/j.envsoft.2021.105069]

11. Binois, M.; Ginsbourger, D.; Roustant, O. Quantifying uncertainty on Pareto fronts with Gaussian process conditional simulations. Eur. J. Oper. Res.; 2015; 243, pp. 386-394. [DOI: https://dx.doi.org/10.1016/j.ejor.2014.07.032]

12. Bassi, M.; de Cursi, E.S.; Pagnacco, E.; Ellaia, R. Statistics of the Pareto front in multi-objective optimization under uncertainties. Lat. Am. J. Solids Struct.; 2018; 15, [DOI: https://dx.doi.org/10.1590/1679-78255018]

13. Kachitvichyanukul, V. Comparison of Three Evolutionary Algorithms: GA, PSO and DE. Ind. Eng. Manag. Syst.; 2012; 11, pp. 215-223. [DOI: https://dx.doi.org/10.7232/iems.2012.11.3.215]

14. Schwaab, M.; Biscaia, E.C.; Monteiro, J.L.; Pinto, J.C. Nonlinear parameter estimation through particle swarm optimization. Chem. Eng. Sci.; 2008; 63, pp. 1542-1552. [DOI: https://dx.doi.org/10.1016/j.ces.2007.11.024]

15. Benyahia, B.; Latifi, M.A.; Fonteix, C.; Pla, F. Emulsion copolymerization of styrene and butyl acrylate in the presence of a chain transfer agent. Part 2: Parameters estimability and confidence regions. Chem. Eng. Sci.; 2013; 90, pp. 110-118. [DOI: https://dx.doi.org/10.1016/j.ces.2012.12.013]

16. Nogueira, I.B.R.; Martins, M.A.F.; Regufe, M.J.; Rodrigues, A.E.; Loureiro, J.M.; Ribeiro, A.M. Big Data-Based Optimization of a Pressure Swing Adsorption Unit for Syngas Purification: On Mapping Uncertainties from a Metaheuristic Technique. Ind. Eng. Chem. Res.; 2020; 59, pp. 14037-14047. [DOI: https://dx.doi.org/10.1021/acs.iecr.0c01155]

17. Chiandussi, G.; Codegone, M.; Ferrero, S.; Varesio, F.E. Comparison of multi-objective optimization methodologies for engineering applications. Comput. Math. Appl.; 2012; 63, pp. 912-942. [DOI: https://dx.doi.org/10.1016/j.camwa.2011.11.057]

18. Srinivas, N.; Deb, K. Multiobjective optimization using nondominated sorting in genetic algorithms. Evol. Comput.; 1994; 148, pp. 148-162.

19. Dipama, J.; Teyssedou, A.; Aubé, F.; Lizon-a-lugrin, L. A grid based multi-objective evolutionary algorithm for the optimization of power plants. Appl. Therm. Eng.; 2010; 30, pp. 807-816. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2009.12.010]

20. Mousa, A.A.; El-shorbagy, M.A.; Abd-el-wahed, W.F. Local search based hybrid particle swarm optimization algorithm for multiobjective optimization Local search based hybrid particle swarm optimization algorithm for multiobjective optimization. Swarm Evol. Comput.; 2018; 3, pp. 1-14. [DOI: https://dx.doi.org/10.1016/j.swevo.2011.11.005]

21. Tanaka, M.; Watanabe, H.; Furukawa, Y.; Tanino, T. GA-based decision support system for multicriteria optimization. Ind. Manag. Eng. Fields; 1995; 2, pp. 1556-1561.

22. Zitzler, E.; Deb, K. Comparison of Multiobjective Evolutionary Algorithms: Empirical Results. Evol. Comput.; 2000; 8, pp. 173-195. [DOI: https://dx.doi.org/10.1162/106365600568202]

23. Coello, C.A.; Lamont, G.B.; Van Veldhuizen, D.A. Evolutionary Algorithms for Solving Multi-Objective Problems; Springer: New York, NY, USA, 2007; ISBN 9780387310299

24. Ranjbar, A.; Talatahari, S.; Hakimpour, F. The application of multi-objective charged system search algorithm for optimization problems. Sci. Iran.; 2019; 26, pp. 12-14. [DOI: https://dx.doi.org/10.24200/sci.2018.20184]

25. Ong, Y.S.; Nair, P.B.; Lum, K.Y. Max-min surrogate-assisted evolutionary algorithm for robust design. IEEE Trans. Evol. Comput.; 2006; 10, pp. 392-404. [DOI: https://dx.doi.org/10.1109/TEVC.2005.859464]

26. Hansen, M.P.; Jaszkiewicz, A. Evaluating the Quality of Approximations to the Non-dominated Set; IMM, Department of Mathematical Modelling, Technical University of Denmark: Lyngby, Denmark, 1998; 31, 30.

27. Deb, K.; Sachin, J. Running Performance Metrics for Evolutionary Multi-Objective Optimization. Kangal Rep.; 2002; 2002004, pp. 13-20.

28. Van Veldhuizen, D.A.; Lamont, G.B. Multiobjective evolutionary algorithms: Analyzing the state-of-the-art. Evol. Comput.; 2000; 8, pp. 125-147. [DOI: https://dx.doi.org/10.1162/106365600568158]

29. Gaspar-Cunha, A.; Antunes, C.H.; Coello, C.C. Evolutionary multi-criterion optimization: 8th international conference, EMO 2015 Guimarães, Portugal, 29 March—1 April 2015 proceedings, Part II. Lect. Notes Comput. Sci.; 2015; 9019, pp. 110-125. [DOI: https://dx.doi.org/10.1007/978-3-319-15892-1]

30. Coello, C.A.C.; Sierra, M.R. A study of the parallelization of a coevolutionary multi-objective evolutionary algorithm. Lect. Notes Artif. Intell.; 2004; 2972, pp. 688-697. [DOI: https://dx.doi.org/10.1007/978-3-540-24694-7_71]

31. Zitzler, E.; Thiele, L. Multiobjective optimization using evolutionary algorithms—A comparative case study. Lect. Notes Comput. Sci.; 1998; 1498, pp. 292-301. [DOI: https://dx.doi.org/10.1007/bfb0056872]

32. Regufe, M.J.; Tamajon, J.; Ribeiro, A.M.; Ferreira, A.F.P.; Lee, U.; Hwang, Y.K.; Chang, J.; Serre, C.; Loureiro, J.M.; Rodrigues, A.E. Syngas Purification by Porous Amino—Functionalized Titanium Terephthalate MIL-125. Energy Fuels; 2015; 29, pp. 4654-4664. [DOI: https://dx.doi.org/10.1021/acs.energyfuels.5b00975]