RELIABILITY-BASED MULTI-OBJECTIVE OPTIMIZATION APPLIED TO CHEMICAL ENGINEERING DESIGN

ABSTRACT

Chemical engineering optimization represents a significant challenge due to the complexity of the mathematical models that are frequently required in this area. These models are normally associated with nonlinear equations that represent mass, energy, and momentum balances, which are submitted to physical, constitutive, environmental, and design limitations. The design of chemical systems is generally carried out by considering the model, the vector of design variables, and system parameters as deterministic values, i.e., small variations in these quantities do not affect the objective function. In this contribution, a new methodology based on a double loop iteration process to evaluate the influence of uncertainties on chemical engineering design is proposed. The inner optimization loop is used to find the solution associated with the highest probability value by using the so-called Inverse Reliability Analysis and the outer loop is the regular optimization loop used to determine the vector of design variables. For this aim, the Multi-Objective Optimization Water Cycle Algorithm is improved, adopting a mechanism of neighborhood exploration. For illustration purposes, the proposed methodology is applied to mathematical functions and to chemical engineering design. The obtained results demonstrate that the proposed strategy represents an interesting alternative to reliability design in chemical engineering.

Key words:
Chemical engineering design; Reliability-based optimization; Water cycle algorithm; Nono- and multi-objective problems; Inverse reliability analysis

INTRODUCTION

Traditionally, in chemical engineering design, the model, the vector of design variables, and the vector of parameters are considered as deterministic quantities, i.e., the possible influence of uncertainties on the resulting design is disregarded. However, small variations in these quantities can affect the vector of objective functions and, consequently, the optimal design. Among the main sources of uncertainties observed in engineering processes, we can cite (Ritto et al., 2008; Leidemer, 2009): i) uncertainties in the parameters of the model, such as geometrical and constitutive parameters; and ii) simplifications adopted in the model formulation. It is important to add that other sources may exist, such as the uncertainty associated with process variability, i.e., disturbances in the variables of the process such as air and/or liquid flows, temperature variations and material compositions. In general, even if the global minimum solution is found, it may be difficult to implement it in practice, e.g., due to the requirement of a high degree of fabrication accuracy. In this context, it is necessary to consider the influence of uncertainties during chemical engineering design to minimize the effect of these quantities on the vector of objective functions.

Various probabilistic methods have been proposed over the years to deal with uncertainties. Among these, reliability-based optimization (RBO) appears as one of the most employed approaches. According to Du and Chen (2004), RBO emphasizes high reliability in the design by ensuring the probabilistic achievement of the considered constraint at a desired level. Traditionally, deterministic analyses around the nominal solution are necessary to evaluate the probability value. Monte Carlo simulation is commonly used for this purpose (Du and Chen, 2004). The major difficultly associated with this approach is the associated computational cost (Carter, 1997; Melchers 1999). To overcome this limitation, different optimization-based approaches propose to determine the probability of failure without using simulation methods (Deb et al., 2009). In these approaches, the vector of random variables is converted from the physical space to the standard space. A different optimization problem is formulated to compute the largest probability of failure, which is simultaneously equated to the desired value (Deb et al., 2009). Thus, the RBO approach can be classified into four main categories based on the approach used to determine the probability of failure (Agarwal, 2004; Deb et al., 2009): i) simulation methods, ii) double-loop methods, iii) decoupled methods, and iv) single-loop methods. According to Aoues and Chateauneuf (2009), the convergence of a given RBO problem depends on both the initial design and the optimization method.

Traditionally, RBO problems are treated in the mono-objective context; then, by defining a reliability parameter, the optimization technique is applied. In this case, the computational cost increases since the optimization procedure is applied many times to evaluate the influence of the reliability parameter on the value of the objective function. In this contribution, a multi-objective problem is defined. The original objective is associated with the maximization of the reliability parameter during the reliability analysis. To solve this multi-objective problem, the Multi-objective Optimization using the Water Cycle Algorithm (e-MOWCA) is proposed. The e-MOWCA algorithm is based on the extension of the Water Cycle Algorithm (WCA) to the multi-objective context through the incorporation of three operators, namely non-dominated sorting, crowding distance, and neighborhood exploration. The influence of uncertainties is evaluated by using Inverse Reliability Analysis (IRA) (Du, 2005). The proposed methodology is based on a double-loop iteration process. In the outer optimization loop, e-MOWCA is carried out to define the vector of design variables, which is necessary to solve the reliability-based multi-objective optimization problem. In the inner optimization loop, the IRA procedure is carried out to find the point with the highest associated probability value, in agreement with the optimum values determined by using the outer optimization loop.

The organization of this paper is the following. The next two sections present reviews on WCA and e-MOWCA, respectively. Subsequently, a description of the RBO and IRA concepts are presented. The proposed methodology is then presented in the section devoted to Methodology. For illustration purposes, the section Numerical Results brings applications involving mathematical functions and three chemical engineering test cases, as follows: i) Heat Exchanger Network Design, ii) Reactor Network Design, and iii) Catalyst Mixing Problem. Finally, the conclusions are outlined in the last section of the present contribution.

WATER CYCLE ALGORITHM (WCA)

WCA, an optimization technique that belongs to the family of bio-inspired optimization algorithms, differs from other evolutionary approaches in the scheme used to generate potential candidates to solve the optimization problem. Proposed by Eskandar et al. (2012), this optimization technique is based on the water cycle process (or hydrological cycle or H₂O cycle) that describes water movement on the earth through three basic mechanisms: evaporation, precipitation, and surface run-off. In general, water moves downhill, creating streams and rivers from the mountains to the sea. Streams and rivers collect water from the rain and other streams on their way downhill. In the rivers and lakes, water evaporates when plants give off water, in accordance with the transpiration process. Clouds are generated from evaporated water taken to the atmosphere. The clouds condense in the colder atmosphere and release the water back in the form of rain, thus creating new streams and rivers. The parameters considered by using the WCA are the following: population size, number of generations, number of rivers and constant evaporation condition. The complete description of the proposed operators to update the population of each generation of the WCA can be found in Eskandar et al. (2012).

The WCA strategy has been successfully tested in various fields of science, such as: engineering system design (Eskandar et al., 2012), optimization of truss structures (Eskandar et al., 2013), determination of optimal number, location, and size of multiple types of distributed generation units in distribution systems (Baghipour et al., 2014), solution of a quadratic assignment problem (Parhizgar and Shiri, 2014), optimization of a photovoltaic energy conversion system under partial shading condition (Sarvi et al., 2014), solution of mathematical test cases considering a new strategy to determine the evaporation rate (Sadollah et al., 2015a), optimization of classical benchmark discrete truss design problems (Sadollah et al., 2015b), solution of differential-algebraic optimal control problems (Carvalho, 2015), solution of non-convex short-term hydrothermal scheduling due to valve point effects (Haroon and Malik, 2016), solution of benchmark functions and the chaos suppression problem using back stepping control (Pahnehkolaei et al. 2017), among others.

MULTI-OBJECTIVE OPTIMIZATION WATER CYCLE ALGORITHM (E-MOWCA)

Due to the success observed in applications in the mono-objective context, the WCA was extended to consider multi-objective optimization problems. Sadollah et al. (2015c,d) proposed the Multi-objective Optimization Water Cycle Algorithm (MOWCA). The strategy of the extended WCA to solve problems with multiple objectives is defined through the incorporation of two operators to the original algorithm, namely the mechanisms of rank ordering and crowding distance. The performance of MOWCA was evaluated by using benchmark functions and a mechanical engineering system design problem (Sadollah et al., 2015c,d).

In the present contribution, a new strategy, e-MOWCA, is proposed to solve multi-objective problems. e-MOWCA incorporates to the algorithm proposed by Sadollah et al. (2015c,d) a new operator, namely, the neighborhood exploration. This operator is dedicated to refining the potential solution candidates in each generation through the mechanism of neighborhood exploration, as proposed by Hu et al. (2005). The e-MOWCA algorithm presents the following structure. An initial population of size N is generated randomly. All dominated solutions are removed from the population through the Fast Non-Dominated Sorting operator (Deb, 2001). The population is sorted into non-dominated fronts F_j (sets of vectors that are non-dominated with respect to each other). This procedure is repeated until each vector becomes a member of a front. Then, three candidates are selected randomly from the population. A child (raindrop) is generated and this process continues until N children are created. Starting from population P₁ of size N, neighbors (Z) are generated for every individual of the population, as follows (Hu et al., 2005):

$$Z(x)=\left(x-D_k\left(g\right)/\mathrm{2,}x+D_k\left(g\right)/2\right)$$ (1)

where D_k (g) (D_k (g)=(k/ξ)(U-L)) is a vector dependent of the generation counter g. ξ is the number of pseudo fronts (curves) defined by the user. The initial maximum neighborhood size for a given population is D_k (0)=[U-L]. L and U represent the lower and upper bounds of the design variables, respectively. The pre-defined number of individuals in each pseudo front is given by (Hu et al., 2005):

$$n_k={\mathrm{rn}}_{k-1},k=\mathrm{2,}\dots ,\xi$$ (2)

where n_k is the number of individuals in the k-th front and r (<1) is the reduction rate. For a given population with N individuals, n_k can be calculated as follows:

$$n_k=N\frac{1-r}{1-r^{R_{\mathrm{pf}}}}{r}^{k-1}$$ (3)

where R_pf is the number of pseudo curves (fronts).

According to Hu et al. (2005) and Lobato and Steffen Jr (2011), if r <1 the first pseudo curve has the highest number of individuals. In this case, each pseudo curve presents an exponentially reducing number of solutions, which emphasizes a local search. Great r values result in more solutions in the last pseudo curve, emphasizing the global search. Consequently, the generated neighbors are classified according to the dominance criterion and only the non-dominated neighbors (P₂) will be associated with P₁ to determine P₃. The population P₃ is then classified according to the dominance criterion. If the number of individuals in the population P₃ is larger than a given value defined by the user, it is truncated according to the Crowding Distance criterion (Deb, 2001). The Crowding Distance describes the density of solutions surrounding a given vector. An infinite Crowding Distance (or an arbitrarily large number for practical purposes) is assigned to the vectors with the smallest or largest values of each objective function. For the remaining vectors, the Crowding Distance is calculated according to:

$${\mathrm{dist}}_{x_i}={\displaystyle \sum }_{j=1}^m\frac{f_{j,i+1}-f_{j,i-1}}{\left|f_{j,\max }-f_{j,\min }\right|}$$ (4)

where f_j corresponds to the j-th objective function and m is the number of objective functions.

In this work, the treatment of constraints is performed by using the Static Penalization Method proposed by Castro (2001). Limited values are attributed to each objective function to play the role of penalization, which guarantees that any non-dominated solution dominates any solution that violates at least one constraint. Any solution that violates a constraint will dominate any solution that presents two constraint violations, and so on. Therefore, layers of solutions are obtained and the number of constraint violations corresponds to the solution rank. For a given constrained problem, the vector containing the objective functions is given by:

$$f\left(x\right)\equiv f\left(x\right)+r_p{n}_{\mathrm{viol}}$$ (5)

where f(x) is the vector of objective functions, r_p is the penalty vector, and n_viol is the number of violated constraints.

Regarding the number of objective function evaluations, e-MOWCA requires more evaluations than the original MOWCA proposed by Sadollah et al. (2015c,d) due to the refinement procedure used during the application of the neighborhood exploration operator (the population is evaluated 2N times in each generation). The adopted refinement procedure accelerates the convergence process, but increases the number of objective function evaluations. In this case, population size and generations number can be adjusted to make e-MOWCA as fast as MOWCA and other multi-objective algorithms.

RELIABILITY-BASED OPTIMIZATION (RBO)

Commonly, the solution of deterministic optimization problems is obtained on a given constraint boundary (or at the intersection of more than one constraint boundary), as presented in Figure 1. In this case, any perturbation in the vector of design variables x₁ and x₂ results in an infeasible solution with a probability distribution around the optimal solution.

Figure 1
Concept of the RBO procedure (Deb et al., 2009).

According to Deb et al. (2009), RBO can be understood as the ability of a system or component to perform its required functions under determined conditions. RBO emphasizes high reliability in the design by ensuring the probabilistic achievement of the considered constraints at a pre-defined level (Du, 2005). The RBO approach in the context of engineering systems is based on the probability of the desired performance of the system to fail (Deb et al., 2009). To find a reliable solution (i.e., a small associated probability to produce an infeasible solution), the true optimal solution must be penalized and an interior solution can be chosen within the feasible region. From a given reliability measure α, it is desired to find a feasible solution that ensures that the probability of having an infeasible solution is given by (1-α). The reliability measure is obtained through the considered uncertainties. Mathematically, the RBO problem can be formulated as (Deb et al., 2009):

$$\min f(x_d,x_r)$$ (6)

subject to:

$$\mathrm{Pr}\left(G_j\left(x_d,x_r\right)\le 0\right)\ge {\alpha }_j,j=\mathrm{1,}\dots ,n_g$$ (7)

$$x_{\mathrm{di}}^l\le x_{\mathrm{di}}\le x_{\mathrm{di}}^u,i=\mathrm{1,}\dots ,n_d$$ (8)

where f and G_j are the objective and constraint functions, respectively, x_d is the vector of design variables (deterministic values associated with the lower x_di^l and upper x_di^u limits), x_r is the vector of random variables, n_g is the number of probabilistic constraints, n_d is the number of design variables, n_r is the number of random variables, and α_j is the desired reliability. Pr is the probability of G_j to be less than or equal to zero (G_j >0 indicates failure). The probability of failure can be defined by a cumulative distribution function, as follows (Deb et al., 2009):

$$\mathrm{Pr}\left(G_j\left(x_d,x_r\right)\le 0\right)={{\displaystyle \int }}^{}\dots {\displaystyle \underset{G_{j\le 0}}{\int }{f}_X\left(x_r\right){\mathrm{dx}}_r}$$ (9)

where f_X is a joint probability density function.

Equation (9) should be evaluated along the design space specified by the vector of inequality constraints to determine the probability of failure Pr. The analytical (or numeric) evaluation of Eq. (9) is expensive due to the specified domain. To overcome this difficulty, the original problem with random variables x_r can be transformed into a similar problem with new random variables (u) by using the so-called Rosenblatt transformation (Rosenblatt,1952). In this approach, considering the new search space (u-space), the most probable point (MPP) for failure is found by locating the minimum distance between the origin and the limit-state (or the constraint function), which is defined by the reliability coefficient β. The probabilistic constraint in Eq. (9) can be further expressed through an inverse transformation, as given by Eq. (10).

$$\mathrm{Pr}\left(G_j\left(x_d,x_r\right)\le 0\right)=\Phi \left(\beta \right)$$ (10)

where Φ is the standard normal distribution function. This transformation is used by different approaches to measure the associated reliability and avoid the evaluation of the integral given by Eq. (9). In the reliability coefficient approach, this parameter is used to describe the probabilistic constraint in Eq. (9).

To solve RBO problems, various approaches are available. Four of the most commonly used are the following: First Order Reliability Method (FORM), Second Order Reliability Method (SORM), Sequential Optimization and Reliability Assessment (SORA), and Inverse Reliability Analysis (IRA) (Du, 2005). In general, the idea is to overcome the computational difficulties by simplifying the integral of f_X (x_r ) in Eq. (9) and determining an approximation for G_j (x_d,x_r ) (Deb et al., 2009).

The IRA approach, a strategy considered as a tool to evaluate the probabilistic constraints, is discussed in the following.

INVERSE RELIABILITY ANALYSIS (IRA)

According to Du (2005), IRA consists in the formulation of an inverse problem for reliability analysis. This can be performed to find the value of G^p , as given by Eq. (11):

$$\mathrm{Pr}\left(G\left(x_d,x_r\right)\le G^p\right)=\alpha$$ (11)

Equation (11) indicates that the probability of the performance function G(x_d ,x_r ) is equal to a given reliability measure α if G(x_d ,x_r ) ≤ G^p , in which G^p is estimated by using FORM. In order to apply the FORM approach, the function G’(x_d ,x_r ) is used as expressed by Eq. (12) and the MPP for Pr(G’_j (x_d,x_r )≤0)= Pr(G’_j (x_d,x_r )≤G^P ) is considered as being u*.

$$G^{\text{'}\left(x_d,x_r\right)}=G\left(x_d,x_r\right)-G^p$$ (12)

From FORM, if the probability α is known, the reliability coefficient is given by:

$$\beta =\left|{\Phi }^{-1}\left(\alpha \right)\right|$$ (13)

where the absolute value is used since the reliability coefficient is the minimum distance between the origin and the limit state (or the constraint function) and is always nonnegative.

Figure 2 shows that the MPP u* is a tangent point of the circle with radius β and the performance function G(x_d,x_r )-G^p =0. Additionally, u* is a point associated with the minimum value of G(x_d,x_r ) on the circle. Therefore, the MPP search for an inverse reliability analysis problem becomes: find the minimum value of G(x_d,x_r ) on the β-circle (or β-sphere for a 3-D problem or β-hyper sphere for a higher dimensional problem).

Figure 2
The Inverse MPP Search (Du, 2005).

The mathematical model for the MPP search is then stated as follows: find the MPP u* where the performance function G(u) is minimized while u* remains on the surface of the β-circle. This procedure can be summarized according the following steps (Du, 2005):

i. Inform the input starting point (k = 0, u₀ = 0, β, and distribution type);

ii. Evaluate G_j (u^k ) and ?G_j (u^k ). Compute a as follows:

$$a^k=\frac{\nabla G_j\left(u^k\right)}{\nabla G_j\left(u^k\right)}$$ (14)

iii. Update u^k+1 = -βa^k and k = k + 1;

iv. Repeat this procedure (steps ii and iii) until convergence (|u^k+1 - u^k | > tolerance) is achieved.

METHODOLOGY

This section presents the e-MOWCA+IRA (solution of multi-objective optimization problems considering reliability) methodology proposed in this paper. As mentioned, the e-MOWCA+IRA approach is based on a double loop iteration process. In the outer loop, the e-MOWCA strategy is carried out as a regular optimization loop to optimize the reliability-based optimization problem. IRA is performed in the inner loop to find the MPP, according to the vector of design variables computed by using e-MOWCA. Figure 3 presents the flowchart regarding the e-MOWCA+IRA strategy.

Figure 3
Flowchart of the e-MOWCA+IRA strategy proposed to solve RBO problems.

The e-MOWCA+IRA strategy is written according to the pseudo-code below:

i. Initially, the RBO problem (in terms of the objective function, constraints, number of variables, reliability index, and types of distribution) and the WCA parameters (population size, number of generations, number of rivers, evaporation condition constant, number of pseudo curves and reduction rate) are defined by the user;

ii. Outer loop: the population of candidates (vector of design variables - x_d ) is generated by using the e-MOWCA strategy. This population consists of individuals (deterministic values) created considering the application of the WCA operators to form a new population, as mentioned earlier. In this case, only the vector of design variables are generated;

iii. Inner loop: for each candidate generated by using the WCA operators, the IRA procedure is evaluated to determine the value of inequality constraints through the determination of the vector of random variables - x_r ;

iv. The values of x_d and x_r are used to evaluate the vector of objective functions and the vector of inequality constraints;

v. The treatment of the constraints is performed through the Penalty Function Method;

vi. After the evaluation of the vector of objective functions, the operators defined in the e-MOWCA strategy are applied. First, the fast non-dominated sorting operator is used to classify the population according to the Pareto dominance. In the following, the neighborhood exploration operator is utilized to explore potential solutions. Finally, the crowding distance operator is applied to reduce the number of non-dominated solutions in each generation;

vii. This process is repeated until convergence is achieved (i.e., the maximum number of generations is found).

In summary, to determine the parameter variations in the proposed methodology, the IRA strategy in association with the user-defined parameters - number of random variables, type of distribution (mean value and variation coefficient) and reliability index - is used to evaluate the probabilistic constraints. After the application of the IRA strategy for each potential candidate, the vector of random variables is obtained and, consequently, the constraints and the vector of objective functions are evaluated. This process continues until all the individuals of the population are considered.

It is worth mentioning that a normal distribution is adopted for the design variables. According to Yin and Chen (2006), this type of distribution is commonly used in engineering applications to describe the linear relationship between the changing variance and the mean value; i.e., σ_x =r_xμ_x , where μ_x and σ_x represent the mean and the standard deviation, respectively, and r_x is the variation coefficient. Rosenblatt (1952) claims that to deal with non-normal distributions, it is necessary to convert the mean and the standard deviation of this type of distribution to an equivalent value with normal distribution. According to Rackwitz-Fiessler’s Two-parameter Equivalent Normal Method (Rosenblatt, 1952), the determination of an equivalent normal distribution at a point of interest X^* can be expressed as follows:

$${\sigma }_x^N=\varphi \left\{{\Phi }^{-1}\left[F_x\left(X^{*}\right)\right]\right\}/f_x\left(X^{*}\right)$$ (15)

$${\mu }_x^N=X^{*}-{\Phi }^{-1}\left[F_x\left(X^{*}\right)\right]{\sigma }_x^N$$ (16)

where ϕ is the Probability Density Function (PDF) of the standard normal distribution. F_x and f_x are the Cumulative Distribution Function (CDF) and PDF of the non-normal distribution of X, respectively. In this case, both the mean value and the standard deviation of a particular type of distribution are transformed into equivalent values with normal distribution.

NUMERICAL RESULTS

To evaluate the ability of the proposed methodology to solve multi-objective optimization problems considering reliability (e-MOWCA), the developed multi-objective algorithm will be initially applied to classical benchmark problems (deterministic test cases). Then, the e-MOWCA strategy will be used to solve chemical engineering systems with different levels of complexity in the presence of uncertainties. It is worth mentioning that both algorithms were implemented by considering the Microsoft Windows^® Platform and Matlab Software^® (2MG RAM). The computational time for each application was evaluated by using a PENTIUM IV (3.2 GHz and 2 GB RAM) microcomputer.

Comparative Study using the e-MOWCA, MOWCA, and NSGA-II Approaches

In the present work, the so-called ZDT functions (Zitzler et al., 2000) are used to test the performance of the proposed multi-objective algorithm, as given by Eq. (17) to Eq. (20).

$$ZDT=\left\{\begin{array}{l}\min f_1\left(x\right)\hfill \\ \min f_2\left(x\right)\equiv g\left(x\right)h\left(f_1\left(x\right),g\left(x\right)\right)\hfill \end{array}\right.$$ (17)

where the corresponding g and h functions are given by:

$$ZD{T}_1=\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ g\left(x_2,x_3,\dots ,x_m\right)=1+9{\displaystyle \sum }_{i=2}^m\frac{x_i}{m-1}\hfill \\ h\left(f_1\left(x\right),g\left(x\right)\right)=1-{\left(f_1(x)/g(x)\right)}^{0.5}\hfill \end{array}\right.$$ (18)

$$ZD{T}_2=\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ g\left(x_2,x_3,\dots ,x_m\right)=1+9{\displaystyle \sum }_{i=2}^m\frac{x_i}{m-1}\hfill \\ h\left(f_1\left(x\right),g\left(x\right)\right)=1-{\left(\frac{f_1\left(x\right)}{g\left(x\right)}\right)}^2\hfill \end{array}\right.$$ (19)

$$ZD{T}_3=\left\{\begin{array}{l}{f}_1\left(x\right)=x_1\hfill \\ g\left(x_2,x_3,\dots ,x_m\right)=1+9{\displaystyle \sum }_{i=2}^m\frac{x_i}{m-1}\hfill \\ h\left(f_1\left(x\right),g\left(x\right)\right)=1-{\left(\frac{f_1\left(x\right)}{g\left(x\right)}\right)}^{0.5}-\left(\frac{f_1\left(x\right)}{g\left(x\right)}\right)\sin \left(10\pi f_1\left(x\right)\right)\hfill \end{array}\right.$$ (20)

The test functions ZDT₁ and ZDT₂ present convex and non-convex Pareto optimal fronts, respectively. ZDT₃ represents the discreteness feature, in which the Pareto optimal front consists of several disjointed continuous convex parts. For all the problems considered, the number of design variables is equal to 15 (m=15 in Eq. (18) to Eq. (20)) and a set of 1000 uniformly spaced Pareto optimal solutions are chosen to compare the convergence (ϒ) and diversity metric (Δ), defined as (Zitzler et al., 2000; Deb, 2001):

$${\rm Y}=\frac{{{\displaystyle \sum }}_{i=1}^{\left|Q\right|}{d}_i}{\left|Q\right|}$$ (21)

$$\Delta =\frac{d_f+d_1+{\displaystyle \sum _{i=1}^{\left|Q\right|-1}}\left|d_i-\overline{d}\right|}{d_f+d_1+\left(\left|Q\right|-1\right)\overline{d}}$$ (22)

where d_i is the Euclidean distance calculated between the Q non-dominated solution obtained by the studied algorithms and the Pareto Optimal solution (analytical solution), is the average distance, d_l and d_f are the distances between the extreme solutions of the Pareto Optimal.

The parameters considered for each multi-objective optimization algorithm are presented in Table 1. The metrics values obtained after 50 runs of each ZDT function are summarized in Table 2. Note that the metrics computed by using the e-MOWCA strategy are similar to those obtained by using NSGA II and MOWCA. In this context, the performance of the proposed methodology considering these metrics is similar that those obtained by other evolutionary algorithms.

Table 1
Parameters considered for the various multi-objective optimization algorithms to solve the ZDT functions.

Parameter	NSGA II (Deb, 2001)	MOWCA (Sadollah et al., 2015c,d)	e-MOWCA
population size	100	100	100
selection strategy	binary tournament	-	-
crossover probability	0.9
mutation probability	0.02
number of generations	200	200	100
number of rivers	-	4	4
evaporation condition constant	-	10-5	10-5
number of pseudo curves	-	-	10
reduction rate	-	-	0.9

Table 2
Average (ϒ, Δ) and variance (σ²) values obtained by the algorithms.

	NSGA II				MOWCA				e-MOWCA
	ϒ	σ2(ϒ)	Δ	σ2(Δ)	ϒ	σ2(ϒ)	Δ	σ2(Δ)	ϒ	σ2(ϒ)	Δ	σ2(Δ)
ZDT1	0.0026	0.0001	0.1540	0.0012	0.0024	0.0001	0.1470	0.0021	0.0022	0.0003	0.1623	0.0033
ZDT2	0.0032	0.0002	0.2130	0.0062	0.0030	0.0001	0.1930	0.0090	0.0029	0.0001	0.1838	0.0041
ZDT3	0.0017	0.0001	0.2351	0.0088	0.0014	0.0002	0.2208	0.0088	0.0012	0.0002	0.2088	0.0013

The number of objective function evaluations in each algorithm is 20100, i.e., all algorithms presented the same computational cost. In terms of computational time, all algorithms required 20 seconds to obtain the Pareto’s Curve, approximately.

Figure 4 presents the Pareto’s Curve for the three ZDT functions considering MOWCA, NSGA II, and e-MOWCA strategies, where it is possible to observe that the results obtained by using e-MOWCA are close to the ones determined by other evolutionary strategies and the analytical solution (Pareto’s Curve).

Figure 4
Solutions obtained considering the three ZDT functions and the different algorithms.

Chemical Engineering Test Cases

In this section, the proposed e-MOWCA+IRA methodology is evaluated considering different chemical engineering design problems: i) Heat Exchanger Network Design, ii) Reactor Network Design, and iii) Catalyst Mixing Problem. It is worth mentioning that originally all the mentioned problems were studied in the deterministic context. In this contribution, these problems will be studied considering uncertainties. For this purpose, the mean values were defined by considering the original mathematical model and the variation coefficients were defined by considering a few runs to evaluate their influence on the optimization process. In the context of the mentioned applications, the following points should be taken into account:

i. The multi-objective problem formulated for each test case consists of the original objective function minimization (or maximization) and the maximization of the reliability coefficient. Thus, each application is formulated as a bi-objective optimization problem;

ii. In order to evaluate the performance of the proposed methodology (e-MOWCA+IRA), the e-MOWCA is associated with another classical strategy to deal with reliability, i.e., the FORM strategy (First Order Reliability Method);

iii. The parameters used by e-MOWCA in both algorithms are presented in Tab. 3. For these parameters, the number of objective function evaluations in e-MOWCA+IRA and e-MOWCA+FORM are 50+100×100×N_IRA and 50+100×100×N_FORM , respectively. N_IRA and N_FORM are the average number of evaluations required by IRA and FORM during the search process, respectively;

Table 3
Parameters considered by using e-MOWCA to solve the chemical engineering test cases.

Parameter	e-MOWCA
population size	50
number of generations	100
number of rivers	4
evaporation condition constant	10-5
number of pseudo curves	10
reduction rate	0.9

iv. The random variables u are considered as being equal to zero in the first iteration of IRA and FORM;

v. The derivatives of G_j were obtained analytically due to their simplicity;

vi. The stopping criterion used to finish the e-MOWCA strategy is the maximum number of generations;

vii. The Euclidean norm of the vector u being smaller than 1.0×10^-5 along two consecutive iterations was considered as the stopping criterion used to conclude IRA and FORM strategies;

viii. For each test case, to formulate the reliability-based multi-objective optimization problem, the components of the vector of random variables (x_r ) are considered as uncertain quantities and the vector of design variables (x_d ) is considered as being composed by deterministic values;

ix. The type of probability distribution, the mean values, and the variation coefficients are defined for each design variable for each test case;

x. For all test cases, the following range for the reliability coefficient (β) is considered: 0.1 ≤ β ≤ 5 (corresponding to the reliability measure (α): 53.9828% ≤ α ≤ 99.9999 %).

Heat Exchanger Network Design

This problem consists in minimizing the overall heat exchange area in a heat exchanger network, as shown in Figure 5 (Floudas and Pardalos, 1990; Angira and Babu, 2003). In this case, the stream should be heated from 100 ^oF to 500 ^oF by using three hot streams with different inlet temperatures.

Figure 5
Heat exchanger network design problem (Angira and Babu, 2003).

Mathematically, this optimization problem can be formulated as follows (Angira and Babu, 2003):

$$\min f=x_1+x_2+x_3$$ (23)

subject to:

$$0.0025(x_4+x_6)-1=0$$ (24)

$$0.0025(-x_4+x_{5+}{x}_7)-1=0$$ (25)

$$0.01(-x_5+x_8)-1=0$$ (26)

$$G_1\equiv 100x_1-x_1{x}_6+833.3325x_4-83333.333\le 0$$ (27)

$$G_2\equiv x_2{x}_4-x_2{x}_7-1250x_4+1250x_5\le 0$$ (28)

$$G_3\equiv x_3{x}_5-x_3{x}_8-2500x_5+1250000\le 0$$ (29)

where x ₁, x ₂, and x ₃ are areas of heat exchangers and x ₄, x ₅, x ₆, x ₇, and x ₈ are temperatures of streams, as shown in Fig. 5. Equations (24) to (26) represent the energy balances for each heat exchanger and Eqs. (27) to (29) represents the maximum demands of heating utilities for each heat exchanger. The design space for this problem is given by (Angira and Babu, 2003): 100≤x ₁≤10000, 1000≤x ₂, x ₃≤10000, and 10≤x ₄, x ₅, x ₆, x ₇, x ₈≤1000.

The global optimum reported by Floudas and Pardalos (1990) and obtained by Angira and Babu (2003) by using the Differential Evolution approach is: [x ₁ x ₂ x ₃ x ₄ x ₅ x ₆ x ₇ x ₈ f] =[579.19 1360.13 5109.92 182.01 295.60 217.92 286.40 395.60 7049.25].

In the reliability context, the multi-objective problem is formulated in terms of the vector of design variables (x _di , i=1, 2, …, 8) (defined as deterministic quantities) and of the vector of random variables (x _rj , j=1, 2, …, 7). Mathematically:

$$\left\{\begin{array}{l}\max \beta \hfill \\ \min f=x_{d1}+x_{d2}+x_{d3}\hfill \end{array}\right.$$ (30)

subject to:

$$0.0025(x_{d4}+x_{d6})-1=0$$ (31)

$$0.0025(-x_{d4}+x_{d5+}{x}_{d7})-1=0$$ (32)

$$0.01(-x_{d5}+x_{d8})-1=0$$ (33)

$$\mathrm{Pr}\left(G_1\le 0\right)\equiv \mathrm{Pr}\left(x_{r1}{x}_{d1}-x_{d1}{x}_{d6}+x_{r2}{x}_{d4}-x_{r3}\le 0\right)=\Phi \left({\beta }_1\right)$$ (34)

$$\mathrm{Pr}\left(G_2\le 0\right)\equiv \mathrm{Pr}\left(x_{d2}{x}_{d4}-x_{d2}{x}_{d7}-x_{r4}{x}_{d4}+x_{r5}{x}_{d5}\le 0\right)=\Phi \left({\beta }_2\right)$$ (35)

$$\mathrm{Pr}\left(G_3\le 0\right)\equiv \mathrm{Pr}\left(x_{d3}{x}_{d5}-x_{d3}{x}_{d8}-x_{r6}{x}_{d5}+x_{r7}\le 0\right)=\Phi \left({\beta }_3\right)$$ (36)

considering the design space as given by: 100≤x _d1 ≤10000, 1000≤x _d2 , x _d3 ≤10000, 10≤x _d4 , x _d5 , x _d6 , x _d7 , x _d8 ≤1000, and 0.1≤β≤5. Reliability is evaluated by considering the mean values and variation coefficients presented in Tab. 4. It is worth mentioning that the mean values were defined by considering the information presented in the original mathematical model. In addition, each variation coefficient was defined after some preliminary runs and its influence on the optimization process.

Figure 6 presents the Pareto’s Curve obtained by using e-MOWCA+IRA and e-MOWCA+FORM considering β₁=β₂=β₃=β in Eqs. (34) to (36) (all constraints present the same reliability coefficient), and normal distribution for all design variables, with the mean values and variation coefficients shown in Tab. 4.

Figure 6
Trade-off frontier between the reliability coefficient and the optimal solution of the heat exchanger network design problem.

Table 4
Statistical data considered in the heat exchanger network design problem.

	Mean value	Variation coefficient		Mean value	Variation coefficient
xr1	100	25 (25%)	xr5	1250	312.5 (25%)
xr2	833.33252	208.33313 (25%)	xr6	2500	625 (25%)
xr3	83333.333	20833.333 (25%)	xr7	1250000	312500 (25%)
xr4	1250	312.5 (25%)

In this figure, note that both algorithms converge to the same Pareto’s Curve, i.e., the IRA strategy was able to obtain a similar result to those obtained by using the FORM strategy. For both algorithms the value of the objective function increases depending on β. Although the design space for the β value has been defined as 0.1≤β≤5, e-MOWCA+IRA and e-MOWCA+FORM converge to values smaller than 1.0. In this case, there are no values of design variables that can satisfy the probabilistic constraints for β>1 (corresponds to 84.1345% of the reliability measure). Consequently, the maximum reliable solution obtained for this test case is 84.1345%. In terms of computational cost, the required average numbers of objective function evaluations were 150050 (N _IRA =15) and 150200 (N _FORM =16), respectively. The deterministic solution results 10050 evaluations. As mentioned, this difference is due to the introduction of IRA (or FORM), which is necessary to evaluate the influence of uncertainties on the optimization process. In terms of computational time, the e-MOWCA+IRA and e-MOWCA+FORM strategies required 150 and 160 seconds to obtain the Pareto’s Curve, respectively.

Table 5 presents some points of Pareto’s Curve (A to F, see Fig. 6) obtained by the proposed methodology (e-MOWCA+IRA). Note that x_d2 and x_d3 increase depending on β (the increase of the reliability coefficient implies an increase of the area of the heat exchanger). Oscillatory behavior is observed for the remaining variables.

Table 5
RBO results associated with the heat exchanger network design problem considering different values of β by using the e-MOWCA+IRA strategy.

	β	f – Eq. (30)	xd1	xd2	xd3	xd4	xd5	xd6	xd7	xd8
	01	7049.266	576.621	1364.201	5108.445	181.793	295.663	218.207	286.130	395.663
A	0.1102	7897.931	605.455	1099.417	6193.058	172.527	271.706	227.473	300.821	371.706
B	0.3182	9275.505	595.334	1555.535	7124.634	165.959	271.533	234.041	294.426	371.533
C	0.5582	10694.517	177.425	1854.990	8662.101	114.303	246.889	285.697	267.414	346.889
D	0.7402	12045.834	290.630	2180.113	9575.090	113.148	244.571	286.852	268.577	344.571
E	0.9122	14793.389	846.742	4009.206	9937.440	143.035	261.849	256.965	281.186	361.849
F	0.9802	22174.691	2429.529	9773.302	9971.858	128.850	274.645	271.150	254.205	374.645

• ¹ Deterministic solution and ² Reliability solution.

The influence of the distribution type on the Pareto’s Curve was also evaluated, in which combinations of the normal (N) and lognormal (LN) distributions were considered as follows: Type I ([x_r1x_r2x_r3x_r4x_r5x_r6x_r7 ]=[N N N N N N N]), Type II ([x_r1x_r2x_r3x_r4x_r5x_r6x_r7 ]=[N N N LN LN LN LN]), and Type III ([x_r1x_r2x_r3x_r4x_r5x_r6x_r7 ]=[LN LN LN N N N N]). Figure 7 presents the obtained Pareto’s Curves by applying e-MOWCA+IRA to solve the heat exchanger network design.

Figure 7
Influence of the type of distribution in trade-off frontier between the reliability coefficient and the optimal solution of the heat exchanger network design problem.

Note that the value of the objective function increases depending on β. Additionally, it can be observed that the results were influenced by the type of distribution considered. Types I, II and III present similar profiles belonging to the interval 0.1≤β≤0.9. In other regions of the specified domain, the e-MOWCA+IRA does not find solutions that satisfy the probabilistic constraints.

Reactor Network Design

Ryoo and Sahinidis (1995) and Angira and Babu (2003) previously studied the next example, which is associated with a reactor network design, as shown in Fig. 8.

Figure 8
Reactor network design problem (Angira and Babu, 2003).

This system represents a sequence of two CSTR (Continuous Stirred Tank Reactor) reactors, where the consecutive reaction A→B→C takes place (A, B and C are the components of this system). The objective of the problem is to maximize the concentration of the product B in the output stream. According to Angira and Babu (2003), this problem becomes complex due to presence of a local minimum solution (f=-0.375), which is close to a global minimum solution (f=-0.3881).

Mathematically, the problem is formulated as follows:

$$\max f=x_4$$ (37)

subject to:

$$x_1+k_1{x}_1{x}_5=1$$ (38)

$$x_2-x_1+k_2{x}_2{x}_6=0$$ (39)

$$x_3+x_1+k_3{x}_3{x}_5=1$$ (40)

$$x_4-x_3+x_2-x_1+k_4{x}_4{x}_6=0$$ (41)

$$G\equiv x_5^{0.5}+x_6^{0.5}-4\le 0$$ (42)

$$\left({\mathrm{0,0,0,0,10}}^{-5}{\mathrm{,10}}^{-5}\right)\le \left(x_1,x_2,x_3,x_4,x_5,x_6\right)\le \left(\mathrm{1,1,1,1,,16,16}\right)$$ (43)

where x _j (j=1, 2, 3, 4) represents the concentration of species A and B in the output stream of reactors 1 and 2, respectively, x _k (k=5, 6) represents the volume of reactors 1 and 2, respectively, and k _i (i=1, …, 4) represents the reaction constants (k ₁=0.09755988, k ₂=0.99k ₁, k ₃=0.0391908 and k ₄=0.9k ₃). Equations (38) to (41) represent the mass balances and Eq. (42) represents an empirical equation that relates the volume of the two reactors as mentioned by Ryoo and Sahinidis (1995).

Angira and Babu (2003) determined the global optimum solution ([x₁x₂x₃x₄x₅x₆f]=[0.771462 0.516997 0.204234 0.388812 3.036504 5.096052 0.388812]) of this problem in the deterministic context by using the Differential Evolution algorithm. In the reliability context, this problem is transformed in a multi-objective problem formulated in terms of the vector of design variables (x_di , i=1, 2, …, 6) and the vector of random variables (x _rj , j=1, 2, …, 4):

$$\left\{\begin{array}{l}\max \beta \hfill \\ \max f=x_{d4}\hfill \end{array}\right.$$ (44)

subject to:

$$x_{d1}+x_{r1}{x}_{d1}{x}_{d5}=1$$ (45)

$$x_{d2}-x_{d1}+x_{r2}{x}_{d2}{x}_{d6}=0$$ (46)

$$x_{d3}+x_{d1}+x_{r3}{x}_{d3}{x}_{d5}=1$$ (47)

$$x_{d4}-x_{d3}+x_{d2}-x_{d1}+x_{r4}{x}_{d4}{x}_{d6}=0$$ (48)

$$\mathrm{Pr}\left(G\le 0\right)\equiv \mathrm{Pr}\left(x_{d5}^{0.5}+x_{d6}^{0.5}-4\le 0\right)=\Phi \left(\beta \right)$$ (49)

$$\left({\mathrm{0,0,0,0,10}}^{-5}{\mathrm{,10}}^{-5}\right)\le \left(x_{d1},x_{d2},x_{d3},x_{d4},x_{d5},x_{d6}\right)\le \left(\mathrm{1,1,1,1,,16,16}\right)$$ (50)

The associated mean values and variation coefficients are shown in Table 6. To analyze the reliability in both algorithms, the mean values were defined considering the parameters presented in the original test case. The variation coefficients were defined after some preliminary runs to evaluate its influence on the optimization process.

Table 6
Statistical data considered in the design of the reactor network.

	Mean value	Variation coefficient
xr1	k1	0.15k1 (15%)
xr2	k2	0.15k2 (15%)
xr3	k3	0.15k3 (15%)
xr4	k4	0.15k4 (15%)

Figure 9 presents the Pareto’s Curve obtained by using the e-MOWCA+IRA and e-MOWCA+FORM strategies considering that all variables follow a normal distribution with the mean values and variation coefficients presented in Table 6. In this figure, it is important to observe that both the algorithms converge to the same Pareto’s Curve, thus evidencing the quality of the proposed strategy to deal with reliability in comparison with the FORM strategy. In addition, the value of the objective function decreases (maximization problem) as β increases. This result is different from the earlier test case, where the design space for β converges for a smaller range as compared with the original design space; namely, the problem of reactor network design converges to values in the range 0.1≤β≤5. In this case, all the values of the design variables were able to satisfy the probabilistic constraint for any value of β. In terms of computational cost, e-MOWCA+IRA required 80050 evaluations of the objective function (mean value; N_IRA =8), e-MOWCA+FORM required 80150 evaluations of the objective function (mean value; N_FORM =9) and 10050 evaluations were necessary for the deterministic context. In terms of computational time, the e-MOWCA+IRA and e-MOWCA+FORM required 80 and 90 seconds to obtain the Pareto’s Curve, respectively.

Figure 9
Trade-off frontier between the reliability coefficient and the optimal solution of the reactor network design.

Table 7 presents selected points of the Pareto’s Curve obtained by e-MOWCA+IRA. Note that x_d2 increases depending on β and other variables exhibit oscillatory behavior.

Table 7
RBO results associated with the reactor network design considering different values of β by using the e-MOWCA+IRA strategy.

	β	f – Eq. (44)	xd1	xd2	xd3	xd4	xd5	xd6
	01	0.38881	0.77152	0.51699	0.20419	0.38881	3.03556	5.09726
A	0.126592	0.38646	0.75384	0.52250	0.21761	0.38646	3.34705	4.58392
B	0.894192	0.37435	0.84755	0.53584	0.14217	0.37435	1.84367	6.02289
C	2.130962	0.35688	0.87178	0.56530	0.12106	0.35688	1.50756	5.61317
D	3.098752	0.34465	0.90876	0.58069	0.08769	0.34465	1.02901	5.84936
E	4.393392	0.32782	0.90502	0.61060	0.09113	0.32782	1.07568	4.99235
F	4.998732	0.31460	0.78559	0.63752	0.19322	0.31460	2.79749	2.40466

• ¹ Deterministic solution and ² Reliability solution.

The influence of the type of distribution on the Pareto’s Curve was also evaluated by considering the following combinations of the normal (N) and lognormal (LN) distributions: Type I ([x _r1 x _r2 x _r3 x _r4 ]=[N N N N]), Type II ([x _r1 x _r2 x _r3 x _r4 ]=[N N LN LN]), Type III ([x _r1 x _r2 x _r3 x _r4 ]=[LN LN N N]), and Type IV ([x _r1 x _r2 x _r3 x _r4 ]=[LN LN LN LN]). Figure 10 presents the Pareto’s Curve as obtained by using e-MOWCA+IRA to solve the problem of reactor network design. For all configurations (Types I, II, III, and IV) the results are similar considering β<1 (approximately). The distributions Type I and II present similar behavior for the remaining values of β. Discrepant behavior is verified for the distributions Type I and II. Therefore, the type of distribution influences the objective function calculation. Additionally, the objective function value decreases as β increases (maximization problem).

Figure 10
Influence of the type of distribution on the trade-off frontier between the reliability coefficient and the optimal solution of the reactor network design.

Catalyst Mixing Problem

The last test case is a classical optimal control problem that considers a plug-flow reactor packed with two catalysts and involving the reactions (Gunn and Thomas, 1965; Logsdon, 1990; Vassiliadis, 1993):

$$S_1\stackrel{k_1,k_2}{\leftrightarrow }{S}_2\stackrel{k_3}{\to }{S}_3$$ (51)

in which the symbols k₁ and k₂ are the reaction rate constants of the first two reactions in a reactor where the catalyst consists entirely of the substance that catalyzes the reversible reactions S₁S₂, while the symbol k₃ is the reaction rate constant of the third reaction S₂→S₃.

The optimal mixing policy of the two catalysts must be determined to maximize the production of species S₃, as follows:

$$\max f=1-X_1(t_f)-X_2(t_f)$$ (52)

subject to mass balance (non-dimensional form) of species S₁ (X₁) and S₂ (X₂):

$$\frac{{\mathrm{dX}}_1}{\mathrm{dt}}=\gamma \left(10X_2-X_1\right)X_1\left(0\right)=1$$ (53)

$$\frac{{\mathrm{dX}}_2}{\mathrm{dt}}=\gamma \left(X_1-10X_2\right)-\left(1-\gamma \right)X_2{X}_2\left(0\right)=0$$ (54)

$$0\le \gamma \le \mathrm{1,}{t}_f=1$$ (55)

where t represents the residence time of the substances from the instant of entry to the reactor and t_f is the final time. The catalyst blending fraction γ (control variable) is the fraction of the catalyst formed by the substance that catalyzes the reaction S₁S₂. This fraction can be varied along the axial position of the reactor.

This classical problem presents the differential index equal to 3, which was first reported by Gunn and Thomas (1965) (f=0.047990) and solved by Logsdon (1990) (f=0.047989) by using orthogonal collocation on finite elements, Vassiliadis (1993) (f=0.048055) by using the control parameterization technique, Lobato et al. (2011) (f=0.048079) by using the Improved Differential Evolution Adaptive with control parameterization technique, etc.

The solution of this problem is performed by discretizing the control variable to transform the original problem with differential index fluctuation into a nonlinear optimization problem with constant differential index (equal to 1). For this purpose, the time interval t which belongs to [t₀t₁t₂ … t_f ] is discretized by using N_elem time nodes. For each subinterval [t_it_i+1 ], i=1, 2, ..., N_elem , the control input is approximated by:

$$\gamma \equiv {\gamma }_i{\mathrm{for\; t}}_i\le t\le t_{i+1}$$ (56)

The unknown control input γ(t) is replaced by N_elem unknown parameters γ₁, γ₂, …, γ_Nelem considering a piecewise linear approximation. In this formulation, the localization of each event t_i is also unknown and should be calculated, resulting in (2N_elem -1) design variables (for this problem, the final time t_f is known). The resulting non-linear optimization problem (with differential index equal to 1) is solved by using the e-MOWCA+IRA and e-MOWCA+FORM strategies.

In the reliability context, the proposed multi-objective problem is formulated in terms of the control variable x_r (all other parameters are defined as deterministic quantities):

$$\left\{\begin{array}{l}\max \beta \hfill \\ \max f=1-X_1(t_f)-X_2(t_f)\hfill \end{array}\right.$$ (57)

subject to:

$$\frac{{\mathrm{dX}}_1}{\mathrm{dt}}=x_r\left(10X_2-X_1\right)X_1\left(0\right)=1$$ (58)

$$\frac{{\mathrm{dX}}_2}{\mathrm{dt}}=x_r\left(X_1-10X_2\right)-\left(1-x_r\right)X_2{X}_2\left(0\right)=0$$ (59)

In this case, the original constraint is evaluated by using the following probabilistic constraint:

$$\mathrm{Pr}\left(G_1\le 0\right)\equiv \mathrm{Pr}\left(-x_r\le 0\right)=\Phi \left({\beta }_1\right)$$ (60)

$$\mathrm{Pr}\left(G_2\le 0\right)\equiv \mathrm{Pr}\left(x_r-1\le 0\right)=\Phi \left({\beta }_2\right)$$ (61)

In addition, the time nodes ([t₀t₁t₂ … t_Nelem-1 ]) are considered deterministic design variables, which are also determined by using the e-MOWCA+IRA and e-MOWCA+FORM strategies:

$$0\le t_i\le \mathrm{1,}i=\mathrm{1,}\dots ,N_{\mathrm{elem}}-1$$ (62)

For this purpose, the associated mean values and variation coefficients for the control variables are presented in Table 8. As explained earlier, the mean values were defined considering the original optimal control problem and the variation coefficients were defined after some preliminary runs to evaluate their influence on the optimization process.

Table 8
Statistical data considered in the catalyst-mixing problem.

	Mean value	Variation coefficients
γ1	xr1	0.1xr1 (10%)
γ2	xr2	0.1xr2 (10%)
⁞	⁞	⁞
γNelem	xrNelem	0.1xrNelem (10%)

Figure 11 presents the Pareto’s Curve obtained by using both strategies considering β₁=β₂=β (all constraints present the same reliability coefficient) and that all variables follow a normal distribution with the mean values and variation coefficients presented in Table 8.

Figure 11
Trade-off frontier between the reliability coefficient and the optimal solution of the catalyst-mixing problem.

In Figure 11, it is important to observe that both strategies were able to obtain the same Pareto’s Curve (demonstrating the quality of the solution obtained by using the IRA and FORM strategies). Note that the value of the objective function decreases (maximization problem) as β increases. In terms of the average number of objective function evaluations, this problem required 120050 evaluations (N _IRA =12), 120050 evaluations (N _FORM =12) and 10050 evaluations for the deterministic case. In terms of computational time, the e-MOWCA+IRA and e-MOWCA+FORM strategies required 250 and 260 seconds to obtain the Pareto’s Curve, respectively.

Table 9 presents selected points of the Pareto’s Curve (A to F, see Figure 11) as obtained by e-MOWCA+IRA. Note that the values of x_r2 and x_r3 increase with β. Differently, x_r1 decreases as β increases to guarantee the obedience of the probabilistic constraints.

Table 9
RBO results associated with the catalyst-mixing problem considering different values of β by using the e-MOWCA+IRA strategy.

	β	f – Eq. (57)	ts1	ts1	xr1	xr2	xr3
	01	0.04807	0.13447	0.73057	1	0.22481	0
A	0.102372	0.04801	0.13748	0.73212	0.99175	0.21209	0.00019
B	1.276152	0.04771	0.15241	0.70588	0.99773	0.28207	0.00198
C	2.677362	0.04723	0.20593	0.74175	0.99992	0.30377	0.00082
D	3.645132	0.04671	0.24998	0.73648	0.99637	0.35480	0.00049
E	4.548212	0.04612	0.27170	0.73719	0.99606	0.41475	0.00177
F	4.999752	0.04567	0.35337	0.70625	0.99551	0.42915	0.00318

• ¹ Deterministic solution and ² Reliability solution.

The influence of the type of distribution on the Pareto’s Curve was evaluated by using the following combinations of the normal (N) and lognormal (LN) distributions: Type I ([x _r1 x _r2 x _r3 ]=[N N N]) and Type II ([x _r1 x _r2 x _r3 ]=[LN LN LN]). The influence of the distribution type by using e-MOWCA+IRA to solve the catalyst mixing problem is presented in Fig. 12. For all configurations, the value of the objective function decreases as β increases. Additionally, both distribution types present the same behavior for β≤1.5 (approximately). However, for β>1.5, the distributions Type I and II exhibit discrepant behavior.

Figure 12
Influence of the type of distribution on the trade-off frontier between the reliability coefficient and the optimal solution of the catalyst-mixing problem.

CONCLUSIONS

In this contribution, a new strategy to solve reliability problems was proposed. The strategy named as e-MOWCA (Multi-objective Optimization Water Cycle Algorithm) is based on the extension of the Water Cycle Algorithm proposed by Sadollah et al. (2015c,d) to the multi-objective context. This extension consists in the incorporation of the neighborhood exploration operator, in association with Inverse Reliability Analysis (IRA) to deal with uncertainties.

The proposed multi-objective approach was applied to solve classical mathematical problems. The obtained results demonstrated that the algorithm was able to obtain the Pareto’s Curve in terms of convergence and diversity. The methodology proposed to deal with uncertainties (e-MOWCA+IRA) was applied to solve chemical engineering design problems that exhibit different complexities. The test cases evaluated in this contribution were solved first by other authors considering the deterministic context only. For this purpose, all the design variables were defined as deterministic quantities, i.e., only the vector of random variables contains uncertainties. In addition, the number of both random and deterministic variables and their corresponding values were defined by the authors as based on each deterministic optimization problem and by considering some preliminary runs to evaluate the influence of these values on the optimization process. Each multi-objective problem was formulated by considering the original objective function associated with the maximization of the reliability coefficient. The results obtained by using the proposed methodology were compared with those obtained by using the e-MOWCA+FORM strategy. In general, similar results (in terms of convergence, number of objective function evaluations and computational time) were obtained by both strategies. This shows the capacity of the IRA strategy to deal with uncertainties in comparison with the FORM strategy.

It is worth mentioning that the type of probability distribution considered influences the solution. As expected, the number of objective function evaluations is greater than the number of evaluations required by classical approaches due to the nature of the adopted optimization strategy.

It is worth mentioning that, from the mathematical point of view, a deterministic solution is more sensitive to small perturbations than a reliable solution, i.e., if both solutions (deterministic and reliable) are perturbed, the mean and variation coefficient will be higher for a deterministic solution as compared to a reliable solution. Therefore, the reliability-based optimized designs are inherently more reliable than the corresponding deterministic solution, which justifies the interest in the methodology presented. Obviously, the procedure requires more objective function evaluations than the deterministic approach.

Although the influence of uncertainties on the vector of design variables was not taken into account, the proposed methodology can be easily extended to more sophisticated design configurations. In addition, a larger number of objective functions can be handled as necessary.

Finally, it can be concluded that the proposed strategy represents an interesting alternative to reliability design of chemical engineering problems.

ACKNOWLEDGEMENT

The authors are thankful to the Brazilian research agencies CNPq, FAPEMIG and CAPES for the financial support of this research work through the National Institute of Science and Technology on Smart Structures for Engineering (INCT-EIE).

REFERENCES

• Agarwal, H., Reliability Based Design Optimization: Formulations and Methodologies. PhD thesis, University of Notre Dame (2004).
• Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Non-linear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 87-91 (2003).
• Aoues, Y., Chateauneuf, A., Benchmark study of numerical methods for reliability-based design. Structural and Multidisciplinary Optimization, 41, 277-294 (2009). https://doi.org/10.1007/s00158-009-0412-2
  » https://doi.org/10.1007/s00158-009-0412-2
• Baghipour, R., Hosseini, S. M., Boor, Z., A Water Cycle Algorithm for Optimal Allocation of DGs in Distribution System Considering Environmental Profit. International Journal of Mechatronics, Electrical and Computer Technology, 4, 430-454 (2014).
• Carvalho, K. T., Solution of Differential-Algebraic Optimal Control Problems with Applications in Engineering. Dissertation, School of Mechanical Engineering, Federal University of Uberlândia (2015).
• Carter, A.D.S., Mechanical Reliability and Design, New York, Wiley (1997). https://doi.org/10.1007/978-1-349-14487-7
  » https://doi.org/10.1007/978-1-349-14487-7
• Castro, R.E., Optimization of Structures with Multi-objective through Genetic Algorithms. Thesis (in Portuguese), COPPE/UFRJ, Brazil (2001).
• Deb, K., .Multi-objective Optimization using Evolutionary Algorithms, John Wiley & Sons, Chichester-UK (2001), ISBN 0-471-87339-X.
• Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliability-based optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 1054-1074 (2009). https://doi.org/10.1109/TEVC.2009.2014361
  » https://doi.org/10.1109/TEVC.2009.2014361
• Du, X., Chen, W., Sequential optimization and reliability assessment method for efficient probabilistic design. Journal of Mechanical Design, 126, 225-233 (2004). https://doi.org/10.1115/1.1649968
  » https://doi.org/10.1115/1.1649968
• Du, X., Probabilistic Engineering Design - First Order and Second Reliability Methods. University of Missouri-Rolla (2005).
• Eskandar, H., Sadollah, A., Bahreininejad, A., Hamdi, M., 2012. Water Cycle Algorithm - A Novel Metaheuristic Optimization Method for Solving Constrained Engineering Optimization Problems. Computers and Structures, 110-111, 151-166 (2012). https://doi.org/10.1016/j.compstruc.2012.07.010
  » https://doi.org/10.1016/j.compstruc.2012.07.010
• Eskandar, H., Sadollah, A, Bahreininejad, A., Weight Optimization of Truss Structures using Water Cycle Algorithm, International Journal of Optimization in Civil Engineering, 3, 115-129 (2013).
• Floudas, C. A., Pardalos, P. M., A Collection of Test Problems for Constrained Global Optimization Algorithms. Lecture notes in computer Science, 455. Springer, Germany (1990). https://doi.org/10.1007/3-540-53032-0
  » https://doi.org/10.1007/3-540-53032-0
• Gunn, D. J., Thomas,W. J., Mass Transport and Chemical Reaction in Multifunctional Catalyst Systems, Chemical Engineering Science, 20, 89-100 (1965). https://doi.org/10.1016/0009-2509(65)85002-3
  » https://doi.org/10.1016/0009-2509(65)85002-3
• Haroon, S. S., Malik, T. N., Evaporation Rate-Based Water Cycle Algorithm for Short-Term Hydrothermal Scheduling, Arabian Journal for Science and Engineering, 1, 1-16 (2016).
• Hu, X., Coello Coello, C.A., Huang, Z.A., New multi-objective evolutionary algorithm: neighborhood exploring evolution strategy, Engineering Optimization, 37, 351-379 (2005). https://doi.org/10.1080/03052150500035658
  » https://doi.org/10.1080/03052150500035658
• Leidemer, M.N., Proposal of Evolutionary Robust Optimization Methodology using the Unscented Transform Applicable to Circuits for RF Circuits/Microwave. Master’s Thesis (in Portuguese), University of Brasília - Brazil (2009).
• Lobato, F.S., Steffen Jr, V., A new multi-objective optimization algorithm based on differential evolution and neighborhood exploring evolution strategy. Journal of Artificial Intelligence and Soft Computing Research, 1, 1-12 (2011).
• Lobato, F.S., Steffen Jr., V., Silva-Neto, A. J., Solution of Singular Optimal Control Problems using the Improved Differential Evolution Algorithm. Journal of Artificial Intelligence and Soft Computing Research , 1, 24-36 (2011).
• Logsdon, J. S., Efficient Determination of Optimal Control Profiles for Differential Algebraic Systems, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA (1990).
• Melchers, R.E., Structural Reliability Analysis and Prediction, John Wiley & Sons, Chichester, England (1999).
• Pahnehkolaei, S. M. A., Alfi, A., Sadollah, A., Kim, J. H., Gradient-based Water Cycle Algorithm with evaporation rate applied to chaos suppression. Applied Soft Computing, 53, 420-440 (2017). https://doi.org/10.1016/j.asoc.2016.12.030
  » https://doi.org/10.1016/j.asoc.2016.12.030
• Parhizgar, M., Shiri, F. M., Solving Quadratic Assignment Problem using Water Cycle Optimization Algorithm. International Journal of Intelligent Information Systems, 3, 75-79 (2014). https://doi.org/10.11648/j.ijiis.s.2014030601.24
  » https://doi.org/10.11648/j.ijiis.s.2014030601.24
• Rosenblatt, M., Remarks on a multivariate transformation. The Annals of Mathematical Statistics, 23, 470-472 (1952). https://doi.org/10.1214/aoms/1177729394
  » https://doi.org/10.1214/aoms/1177729394
• Ritto, T.G., Sampaio, R., Cataldo, E., Timoshenko beam with uncertainty on the boundary conditions. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 30(4), 295-303 (2008). https://doi.org/10.1590/S1678-58782008000400005
  » https://doi.org/10.1590/S1678-58782008000400005
• Ryoo, H. S., Sahinidis, B. P., Global Optimization of Nonconvex NLPs and MINLPs with Application in Process Design. Computers & Chemical Engineering, 19(5), 551-566 (1995). https://doi.org/10.1016/0098-1354(94)00097-2
  » https://doi.org/10.1016/0098-1354(94)00097-2
• Sadollah, A., Eskandar, H., Bahreininejad, A., Kim, J. H., Water cycle algorithm with evaporation rate for solving constrained and unconstrained optimization problems, Applied Soft Computing , 30, 58-71 (2015a). https://doi.org/10.1016/j.asoc.2015.01.050
  » https://doi.org/10.1016/j.asoc.2015.01.050
• Sadollah, A., Eskandar, H., Bahreininejad, A., Kim, J. H., Water cycle, mine blast and improved mine blast algorithms for discrete sizing optimization of truss structures Computers and Structures, 149, 1-16 (2015b). https://doi.org/10.1016/j.compstruc.2014.12.003
  » https://doi.org/10.1016/j.compstruc.2014.12.003
• Sadollah, A., Eskandar, H., Kim, J. H., Bahreininejad, A., Water cycle algorithm for solving multi-objective optimization problems, Soft Computing, 19(9), 2587-2603 (2015c). https://doi.org/10.1007/s00500-014-1424-4
  » https://doi.org/10.1007/s00500-014-1424-4
• Sadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multi-objective problems, Applied Soft Computing , 27, 279-298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
  » https://doi.org/10.1016/j.asoc.2014.10.042
• Sarvi, M., Soltani, I., Avanaki, I. N., A Water Cycle Algorithm Maximum Power Point Tracker for Photovoltaic energy conversion system under Partial Shading Condition. Applied Mathematics in Engineering, Management and Technology, 2, 103-116 (2014).
• Vassiliadis, V., Computational Solution of Dynamic Optimization Problems with General Differential-Algebraic Constraints, Ph.D. thesis, University of London, London, UK (1993).
• Yin, X., Chen, W., Enhanced sequential optimization and reliability assessment method for probabilistic optimization with varying design variance. Structure and Infrastructure Engineering, 2(3-4), 261-275 (2006). https://doi.org/10.1080/15732470600590317
  » https://doi.org/10.1080/15732470600590317
• Zitzler, E., Deb, K., Thiele, L., Comparison of Multiobjective Evolutionary Algorithms: Empirical Results, Evolutionary Computation Journal, 8, 125-148 (2000). https://doi.org/10.1162/106365600568202
  » https://doi.org/10.1162/106365600568202

Publication Dates

• Publication in this collection
  15 July 2019
• Date of issue
  Jan-Mar 2019

History

• Received
  21 July 2017
• Reviewed
  05 Mar 2018
• Accepted
  10 Mar 2018