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 socalled Inverse Reliability Analysis and the outer loop is the regular optimization loop used to determine the vector of design variables. For this aim, the MultiObjective 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; Reliabilitybased optimization; Water cycle algorithm; Nono and multiobjective 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., 2008Ritto, 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), 295303 (2008). https://doi.org/10.1590/S167858782008000400005
https://doi.org/10.1590/S16785878200800...
; Leidemer, 2009Leidemer, 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).): 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, reliabilitybased optimization (RBO) appears as one of the most employed approaches. According to Du and Chen (2004Du, X., Chen, W., Sequential optimization and reliability assessment method for efficient probabilistic design. Journal of Mechanical Design, 126, 225233 (2004). https://doi.org/10.1115/1.1649968
https://doi.org/10.1115/1.1649968...
), 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, 1997Carter, A.D.S., Mechanical Reliability and Design, New York, Wiley (1997). https://doi.org/10.1007/9781349144877
https://doi.org/10.1007/978134914487...
; Melchers 1999Melchers, R.E., Structural Reliability Analysis and Prediction, John Wiley & Sons, Chichester, England (1999).). To overcome this limitation, different optimizationbased approaches propose to determine the probability of failure without using simulation methods (Deb et al., 2009Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliabilitybased optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 10541074 (2009). https://doi.org/10.1109/TEVC.2009.2014361
https://doi.org/10.1109/TEVC.2009.201436...
). 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, 2004Agarwal, H., Reliability Based Design Optimization: Formulations and Methodologies. PhD thesis, University of Notre Dame (2004).; Deb et al., 2009): i) simulation methods, ii) doubleloop methods, iii) decoupled methods, and iv) singleloop methods. According to Aoues and Chateauneuf (2009Aoues, Y., Chateauneuf, A., Benchmark study of numerical methods for reliabilitybased design. Structural and Multidisciplinary Optimization, 41, 277294 (2009). https://doi.org/10.1007/s0015800904122
https://doi.org/10.1007/s001580090412...
), the convergence of a given RBO problem depends on both the initial design and the optimization method.
Traditionally, RBO problems are treated in the monoobjective 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 multiobjective problem is defined. The original objective is associated with the maximization of the reliability parameter during the reliability analysis. To solve this multiobjective problem, the Multiobjective Optimization using the Water Cycle Algorithm (eMOWCA) is proposed. The eMOWCA algorithm is based on the extension of the Water Cycle Algorithm (WCA) to the multiobjective context through the incorporation of three operators, namely nondominated sorting, crowding distance, and neighborhood exploration. The influence of uncertainties is evaluated by using Inverse Reliability Analysis (IRA) (Du, 2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (2005).). The proposed methodology is based on a doubleloop iteration process. In the outer optimization loop, eMOWCA is carried out to define the vector of design variables, which is necessary to solve the reliabilitybased multiobjective 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 eMOWCA, 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 bioinspired 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. (2012Eskandar, 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, 110111, 151166 (2012). https://doi.org/10.1016/j.compstruc.2012.07.010
https://doi.org/10.1016/j.compstruc.2012...
), this optimization technique is based on the water cycle process (or hydrological cycle or H_{2}O cycle) that describes water movement on the earth through three basic mechanisms: evaporation, precipitation, and surface runoff. 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., 2012Eskandar, 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, 110111, 151166 (2012). https://doi.org/10.1016/j.compstruc.2012.07.010
https://doi.org/10.1016/j.compstruc.2012...
), optimization of truss structures (Eskandar et al., 2013Eskandar, H., Sadollah, A, Bahreininejad, A., Weight Optimization of Truss Structures using Water Cycle Algorithm, International Journal of Optimization in Civil Engineering, 3, 115129 (2013).), determination of optimal number, location, and size of multiple types of distributed generation units in distribution systems (Baghipour et al., 2014Baghipour, 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, 430454 (2014).), solution of a quadratic assignment problem (Parhizgar and Shiri, 2014Parhizgar, M., Shiri, F. M., Solving Quadratic Assignment Problem using Water Cycle Optimization Algorithm. International Journal of Intelligent Information Systems, 3, 7579 (2014). https://doi.org/10.11648/j.ijiis.s.2014030601.24
https://doi.org/10.11648/j.ijiis.s.20140...
), optimization of a photovoltaic energy conversion system under partial shading condition (Sarvi et al., 2014Sarvi, 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, 103116 (2014).), solution of mathematical test cases considering a new strategy to determine the evaporation rate (Sadollah et al., 2015Sadollah, 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, 5871 (2015a). https://doi.org/10.1016/j.asoc.2015.01.050
https://doi.org/10.1016/j.asoc.2015.01.0...
a), optimization of classical benchmark discrete truss design problems (Sadollah et al., 2015bSadollah, 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, 116 (2015b). https://doi.org/10.1016/j.compstruc.2014.12.003
https://doi.org/10.1016/j.compstruc.2014...
), solution of differentialalgebraic optimal control problems (Carvalho, 2015Carvalho, K. T., Solution of DifferentialAlgebraic Optimal Control Problems with Applications in Engineering. Dissertation, School of Mechanical Engineering, Federal University of Uberlândia (2015).), solution of nonconvex shortterm hydrothermal scheduling due to valve point effects (Haroon and Malik, 2016Haroon, S. S., Malik, T. N., Evaporation RateBased Water Cycle Algorithm for ShortTerm Hydrothermal Scheduling, Arabian Journal for Science and Engineering, 1, 116 (2016).), solution of benchmark functions and the chaos suppression problem using back stepping control (Pahnehkolaei et al. 2017Pahnehkolaei, S. M. A., Alfi, A., Sadollah, A., Kim, J. H., Gradientbased Water Cycle Algorithm with evaporation rate applied to chaos suppression. Applied Soft Computing, 53, 420440 (2017). https://doi.org/10.1016/j.asoc.2016.12.030
https://doi.org/10.1016/j.asoc.2016.12.0...
), among others.
MULTIOBJECTIVE OPTIMIZATION WATER CYCLE ALGORITHM (EMOWCA)
Due to the success observed in applications in the monoobjective context, the WCA was extended to consider multiobjective optimization problems. Sadollah et al. (2015Sadollah, A., Eskandar, H., Kim, J. H., Bahreininejad, A., Water cycle algorithm for solving multiobjective optimization problems, Soft Computing, 19(9), 25872603 (2015c). https://doi.org/10.1007/s0050001414244
https://doi.org/10.1007/s005000141424...
c,dSadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
https://doi.org/10.1016/j.asoc.2014.10.0...
) proposed the Multiobjective 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., 2015cSadollah, A., Eskandar, H., Kim, J. H., Bahreininejad, A., Water cycle algorithm for solving multiobjective optimization problems, Soft Computing, 19(9), 25872603 (2015c). https://doi.org/10.1007/s0050001414244
https://doi.org/10.1007/s005000141424...
,dSadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
https://doi.org/10.1016/j.asoc.2014.10.0...
).
In the present contribution, a new strategy, eMOWCA, is proposed to solve multiobjective problems. eMOWCA incorporates to the algorithm proposed by Sadollah et al. (2015Sadollah, A., Eskandar, H., Kim, J. H., Bahreininejad, A., Water cycle algorithm for solving multiobjective optimization problems, Soft Computing, 19(9), 25872603 (2015c). https://doi.org/10.1007/s0050001414244
https://doi.org/10.1007/s005000141424...
c,dSadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
https://doi.org/10.1016/j.asoc.2014.10.0...
) 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. (2005Hu, X., Coello Coello, C.A., Huang, Z.A., New multiobjective evolutionary algorithm: neighborhood exploring evolution strategy, Engineering Optimization, 37, 351379 (2005). https://doi.org/10.1080/03052150500035658
https://doi.org/10.1080/0305215050003565...
). The eMOWCA 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 NonDominated Sorting operator (Deb, 2001Deb, K., .Multiobjective Optimization using Evolutionary Algorithms, John Wiley & Sons, ChichesterUK (2001), ISBN 047187339X.). The population is sorted into nondominated fronts F_{j} (sets of vectors that are nondominated 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_{1} of size N, neighbors (Z) are generated for every individual of the population, as follows (Hu et al., 2005):
where D_{k} (g) (D_{k} (g)=(k/ξ)(UL)) 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)=[UL]. L and U represent the lower and upper bounds of the design variables, respectively. The predefined number of individuals in each pseudo front is given by (Hu et al., 2005Hu, X., Coello Coello, C.A., Huang, Z.A., New multiobjective evolutionary algorithm: neighborhood exploring evolution strategy, Engineering Optimization, 37, 351379 (2005). https://doi.org/10.1080/03052150500035658
https://doi.org/10.1080/0305215050003565...
):
where n_{k} is the number of individuals in the kth front and r (<1) is the reduction rate. For a given population with N individuals, n_{k} can be calculated as follows:
where R_{pf} is the number of pseudo curves (fronts).
According to Hu et al. (2005Hu, X., Coello Coello, C.A., Huang, Z.A., New multiobjective evolutionary algorithm: neighborhood exploring evolution strategy, Engineering Optimization, 37, 351379 (2005). https://doi.org/10.1080/03052150500035658
https://doi.org/10.1080/0305215050003565...
) and Lobato and Steffen Jr (2011Lobato, F.S., Steffen Jr., V., SilvaNeto, A. J., Solution of Singular Optimal Control Problems using the Improved Differential Evolution Algorithm. Journal of Artificial Intelligence and Soft Computing Research , 1, 2436 (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 nondominated neighbors (P_{2}) will be associated with P_{1} to determine P_{3}. The population P_{3} is then classified according to the dominance criterion. If the number of individuals in the population P_{3} is larger than a given value defined by the user, it is truncated according to the Crowding Distance criterion (Deb, 2001Deb, K., .Multiobjective Optimization using Evolutionary Algorithms, John Wiley & Sons, ChichesterUK (2001), ISBN 047187339X.). 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:
where f_{j} corresponds to the jth 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 (2001Castro, R.E., Optimization of Structures with Multiobjective through Genetic Algorithms. Thesis (in Portuguese), COPPE/UFRJ, Brazil (2001).). Limited values are attributed to each objective function to play the role of penalization, which guarantees that any nondominated 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:
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, eMOWCA requires more evaluations than the original MOWCA proposed by Sadollah et al. (2015Sadollah, A., Eskandar, H., Kim, J. H., Bahreininejad, A., Water cycle algorithm for solving multiobjective optimization problems, Soft Computing, 19(9), 25872603 (2015c). https://doi.org/10.1007/s0050001414244
https://doi.org/10.1007/s005000141424...
c,dSadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
https://doi.org/10.1016/j.asoc.2014.10.0...
) 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 eMOWCA as fast as MOWCA and other multiobjective algorithms.
RELIABILITYBASED 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_{1} and x_{2} results in an infeasible solution with a probability distribution around the optimal solution.
According to Deb et al. (2009Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliabilitybased optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 10541074 (2009). https://doi.org/10.1109/TEVC.2009.2014361
https://doi.org/10.1109/TEVC.2009.201436...
), 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 predefined level (Du, 2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (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):
subject to:
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., 2009Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliabilitybased optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 10541074 (2009). https://doi.org/10.1109/TEVC.2009.2014361
https://doi.org/10.1109/TEVC.2009.201436...
):
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 socalled Rosenblatt transformation (Rosenblatt,1952Rosenblatt, M., Remarks on a multivariate transformation. The Annals of Mathematical Statistics, 23, 470472 (1952). https://doi.org/10.1214/aoms/1177729394
https://doi.org/10.1214/aoms/1177729394...
). In this approach, considering the new search space (uspace), the most probable point (MPP) for failure is found by locating the minimum distance between the origin and the limitstate (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).
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, 2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (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., 2009Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliabilitybased optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 10541074 (2009). https://doi.org/10.1109/TEVC.2009.2014361
https://doi.org/10.1109/TEVC.2009.201436...
).
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 (2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (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):
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*.
From FORM, if the probability α is known, the reliability coefficient is given by:
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 3D problem or βhyper sphere for a higher dimensional problem).
The Inverse MPP Search (Du, 2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (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, 2005Du, X., Probabilistic Engineering Design  First Order and Second Reliability Methods. University of MissouriRolla (2005).):
i. Inform the input starting point (k = 0, u_{0} = 0, β, and distribution type);
ii. Evaluate G_{j} (u^{k} ) and ?G_{j} (u^{k} ). Compute a as follows:
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 eMOWCA+IRA (solution of multiobjective optimization problems considering reliability) methodology proposed in this paper. As mentioned, the eMOWCA+IRA approach is based on a double loop iteration process. In the outer loop, the eMOWCA strategy is carried out as a regular optimization loop to optimize the reliabilitybased optimization problem. IRA is performed in the inner loop to find the MPP, according to the vector of design variables computed by using eMOWCA. Figure 3 presents the flowchart regarding the eMOWCA+IRA strategy.
The eMOWCA+IRA strategy is written according to the pseudocode 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 eMOWCA 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 eMOWCA strategy are applied. First, the fast nondominated 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 nondominated 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 userdefined 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 (2006Yin, X., Chen, W., Enhanced sequential optimization and reliability assessment method for probabilistic optimization with varying design variance. Structure and Infrastructure Engineering, 2(34), 261275 (2006). https://doi.org/10.1080/15732470600590317
https://doi.org/10.1080/1573247060059031...
), 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 (1952Rosenblatt, M., Remarks on a multivariate transformation. The Annals of Mathematical Statistics, 23, 470472 (1952). https://doi.org/10.1214/aoms/1177729394
https://doi.org/10.1214/aoms/1177729394...
) claims that to deal with nonnormal 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 RackwitzFiessler’s Twoparameter Equivalent Normal Method (Rosenblatt, 1952), the determination of an equivalent normal distribution at a point of interest X^{*} can be expressed as follows:
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 nonnormal 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 multiobjective optimization problems considering reliability (eMOWCA), the developed multiobjective algorithm will be initially applied to classical benchmark problems (deterministic test cases). Then, the eMOWCA 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 eMOWCA, MOWCA, and NSGAII Approaches
In the present work, the socalled ZDT functions (Zitzler et al., 2000Zitzler, E., Deb, K., Thiele, L., Comparison of Multiobjective Evolutionary Algorithms: Empirical Results, Evolutionary Computation Journal, 8, 125148 (2000). https://doi.org/10.1162/106365600568202
https://doi.org/10.1162/106365600568202...
) are used to test the performance of the proposed multiobjective algorithm, as given by Eq. (17) to Eq. (20).
where the corresponding g and h functions are given by:
The test functions ZDT_{1} and ZDT_{2} present convex and nonconvex Pareto optimal fronts, respectively. ZDT_{3} 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., 2000Zitzler, E., Deb, K., Thiele, L., Comparison of Multiobjective Evolutionary Algorithms: Empirical Results, Evolutionary Computation Journal, 8, 125148 (2000). https://doi.org/10.1162/106365600568202
https://doi.org/10.1162/106365600568202...
; Deb, 2001Deb, K., .Multiobjective Optimization using Evolutionary Algorithms, John Wiley & Sons, ChichesterUK (2001), ISBN 047187339X.):
where d_{i} is the Euclidean distance calculated between the Q nondominated 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 multiobjective 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 eMOWCA 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.
Parameters considered for the various multiobjective optimization algorithms to solve the ZDT functions.
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 eMOWCA strategies, where it is possible to observe that the results obtained by using eMOWCA are close to the ones determined by other evolutionary strategies and the analytical solution (Pareto’s Curve).
Chemical Engineering Test Cases
In this section, the proposed eMOWCA+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 multiobjective 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 biobjective optimization problem;
ii. In order to evaluate the performance of the proposed methodology (eMOWCA+IRA), the eMOWCA is associated with another classical strategy to deal with reliability, i.e., the FORM strategy (First Order Reliability Method);
iii. The parameters used by eMOWCA in both algorithms are presented in Tab. 3. For these parameters, the number of objective function evaluations in eMOWCA+IRA and eMOWCA+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;
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 eMOWCA 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 reliabilitybased multiobjective 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, 1990Floudas, 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/3540530320
https://doi.org/10.1007/3540530320...
; Angira and Babu, 2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).). In this case, the stream should be heated from 100 ^{o}F to 500 ^{o}F by using three hot streams with different inlet temperatures.
Heat exchanger network design problem (Angira and Babu, 2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).).
Mathematically, this optimization problem can be formulated as follows (Angira and Babu, 2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).):
subject to:
where x _{1}, x _{2}, and x _{3} are areas of heat exchangers and x _{4}, x _{5}, x _{6}, x _{7}, and x _{8} 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, 2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).): 100≤x _{1}≤10000, 1000≤x _{2}, x _{3}≤10000, and 10≤x _{4}, x _{5}, x _{6}, x _{7}, x _{8}≤1000.
The global optimum reported by Floudas and Pardalos (1990Floudas, 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/3540530320
https://doi.org/10.1007/3540530320...
) and obtained by Angira and Babu (2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).) by using the Differential Evolution approach is: [x
_{1}
x
_{2}
x
_{3}
x
_{4}
x
_{5}
x
_{6}
x
_{7}
x
_{8}
f] =[579.19 1360.13 5109.92 182.01 295.60 217.92 286.40 395.60 7049.25].
In the reliability context, the multiobjective 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:
subject to:
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 eMOWCA+IRA and eMOWCA+FORM considering β_{1}=β_{2}=β_{3}=β 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.
Tradeoff frontier between the reliability coefficient and the optimal solution of the heat exchanger network design problem.
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, eMOWCA+IRA and eMOWCA+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 eMOWCA+IRA and eMOWCA+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 (eMOWCA+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.
RBO results associated with the heat exchanger network design problem considering different values of β by using the eMOWCA+IRA strategy.
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_{r1}x_{r2}x_{r3}x_{r4}x_{r5}x_{r6}x_{r7} ]=[N N N N N N N]), Type II ([x_{r1}x_{r2}x_{r3}x_{r4}x_{r5}x_{r6}x_{r7} ]=[N N N LN LN LN LN]), and Type III ([x_{r1}x_{r2}x_{r3}x_{r4}x_{r5}x_{r6}x_{r7} ]=[LN LN LN N N N N]). Figure 7 presents the obtained Pareto’s Curves by applying eMOWCA+IRA to solve the heat exchanger network design.
Influence of the type of distribution in tradeoff 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 eMOWCA+IRA does not find solutions that satisfy the probabilistic constraints.
Reactor Network Design
Ryoo and Sahinidis (1995Ryoo, H. S., Sahinidis, B. P., Global Optimization of Nonconvex NLPs and MINLPs with Application in Process Design. Computers & Chemical Engineering, 19(5), 551566 (1995). https://doi.org/10.1016/00981354(94)000972
https://doi.org/10.1016/00981354(94)000...
) and Angira and Babu (2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).) previously studied the next example, which is associated with a reactor network design, as shown in Fig. 8.
Reactor network design problem (Angira and Babu, 2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (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 (2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (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:
subject to:
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
_{1}=0.09755988, k
_{2}=0.99k
_{1}, k
_{3}=0.0391908 and k
_{4}=0.9k
_{3}). 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 (1995Ryoo, H. S., Sahinidis, B. P., Global Optimization of Nonconvex NLPs and MINLPs with Application in Process Design. Computers & Chemical Engineering, 19(5), 551566 (1995). https://doi.org/10.1016/00981354(94)000972
https://doi.org/10.1016/00981354(94)000...
).
Angira and Babu (2003Angira, R., Babu, B. V., Evolutionary Computation for Global Optimization of Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).) determined the global optimum solution ([x_{1}x_{2}x_{3}x_{4}x_{5}x_{6}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 multiobjective 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):
subject to:
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.
Figure 9 presents the Pareto’s Curve obtained by using the eMOWCA+IRA and eMOWCA+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, eMOWCA+IRA required 80050 evaluations of the objective function (mean value; N_{IRA} =8), eMOWCA+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 eMOWCA+IRA and eMOWCA+FORM required 80 and 90 seconds to obtain the Pareto’s Curve, respectively.
Tradeoff 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 eMOWCA+IRA. Note that x_{d2} increases depending on β and other variables exhibit oscillatory behavior.
RBO results associated with the reactor network design considering different values of β by using the eMOWCA+IRA strategy.
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 eMOWCA+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).
Influence of the type of distribution on the tradeoff 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 plugflow reactor packed with two catalysts and involving the reactions (Gunn and Thomas, 1965Gunn, D. J., Thomas,W. J., Mass Transport and Chemical Reaction in Multifunctional Catalyst Systems, Chemical Engineering Science, 20, 89100 (1965). https://doi.org/10.1016/00092509(65)850023
https://doi.org/10.1016/00092509(65)850...
; Logsdon, 1990Logsdon, J. S., Efficient Determination of Optimal Control Profiles for Differential Algebraic Systems, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA (1990).; Vassiliadis, 1993Vassiliadis, V., Computational Solution of Dynamic Optimization Problems with General DifferentialAlgebraic Constraints, Ph.D. thesis, University of London, London, UK (1993).):
in which the symbols k_{1} and k_{2} 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_{1}S_{2}, while the symbol k_{3} is the reaction rate constant of the third reaction S_{2}→S_{3}.
The optimal mixing policy of the two catalysts must be determined to maximize the production of species S_{3}, as follows:
subject to mass balance (nondimensional form) of species S_{1} (X_{1}) and S_{2} (X_{2}):
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_{1}S_{2}. 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 (1965Gunn, D. J., Thomas,W. J., Mass Transport and Chemical Reaction in Multifunctional Catalyst Systems, Chemical Engineering Science, 20, 89100 (1965). https://doi.org/10.1016/00092509(65)850023
https://doi.org/10.1016/00092509(65)850...
) (f=0.047990) and solved by Logsdon (1990Logsdon, J. S., Efficient Determination of Optimal Control Profiles for Differential Algebraic Systems, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA (1990).) (f=0.047989) by using orthogonal collocation on finite elements, Vassiliadis (1993Vassiliadis, V., Computational Solution of Dynamic Optimization Problems with General DifferentialAlgebraic Constraints, Ph.D. thesis, University of London, London, UK (1993).) (f=0.048055) by using the control parameterization technique, Lobato et al. (2011Lobato, F.S., Steffen Jr., V., SilvaNeto, A. J., Solution of Singular Optimal Control Problems using the Improved Differential Evolution Algorithm. Journal of Artificial Intelligence and Soft Computing Research , 1, 2436 (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_{0}t_{1}t_{2} … t_{f} ] is discretized by using N_{elem} time nodes. For each subinterval [t_{i}t_{i+1} ], i=1, 2, ..., N_{elem} , the control input is approximated by:
The unknown control input γ(t) is replaced by N_{elem} unknown parameters γ_{1}, γ_{2}, …, γ_{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 nonlinear optimization problem (with differential index equal to 1) is solved by using the eMOWCA+IRA and eMOWCA+FORM strategies.
In the reliability context, the proposed multiobjective problem is formulated in terms of the control variable x_{r} (all other parameters are defined as deterministic quantities):
subject to:
In this case, the original constraint is evaluated by using the following probabilistic constraint:
In addition, the time nodes ([t_{0}t_{1}t_{2} … t_{Nelem1} ]) are considered deterministic design variables, which are also determined by using the eMOWCA+IRA and eMOWCA+FORM strategies:
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.
Figure 11 presents the Pareto’s Curve obtained by using both strategies considering β_{1}=β_{2}=β (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.
Tradeoff frontier between the reliability coefficient and the optimal solution of the catalystmixing 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 eMOWCA+IRA and eMOWCA+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 eMOWCA+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.
RBO results associated with the catalystmixing problem considering different values of β by using the eMOWCA+IRA strategy.
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 eMOWCA+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.
Influence of the type of distribution on the tradeoff frontier between the reliability coefficient and the optimal solution of the catalystmixing problem.
CONCLUSIONS
In this contribution, a new strategy to solve reliability problems was proposed. The strategy named as eMOWCA (Multiobjective Optimization Water Cycle Algorithm) is based on the extension of the Water Cycle Algorithm proposed by Sadollah et al. (2015Sadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (2015d). https://doi.org/10.1016/j.asoc.2014.10.042
https://doi.org/10.1016/j.asoc.2014.10.0...
c,d) to the multiobjective 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 multiobjective 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 (eMOWCA+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 multiobjective 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 eMOWCA+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 reliabilitybased 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 (INCTEIE).
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 Nonlinear Chemical Engineering Process. In: Proceedings of International Symposium on Process Systems Engineering and Control, Mumbai, 8791 (2003).
 Aoues, Y., Chateauneuf, A., Benchmark study of numerical methods for reliabilitybased design. Structural and Multidisciplinary Optimization, 41, 277294 (2009). https://doi.org/10.1007/s0015800904122
» https://doi.org/10.1007/s0015800904122  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, 430454 (2014).
 Carvalho, K. T., Solution of DifferentialAlgebraic 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/9781349144877
» https://doi.org/10.1007/9781349144877  Castro, R.E., Optimization of Structures with Multiobjective through Genetic Algorithms. Thesis (in Portuguese), COPPE/UFRJ, Brazil (2001).
 Deb, K., .Multiobjective Optimization using Evolutionary Algorithms, John Wiley & Sons, ChichesterUK (2001), ISBN 047187339X.
 Deb, K., Padmanabhan, D., Gupta, S., Mall, A. K., Handling uncertainties through reliabilitybased optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13(5), 10541074 (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, 225233 (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 MissouriRolla (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, 110111, 151166 (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, 115129 (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/3540530320
» https://doi.org/10.1007/3540530320  Gunn, D. J., Thomas,W. J., Mass Transport and Chemical Reaction in Multifunctional Catalyst Systems, Chemical Engineering Science, 20, 89100 (1965). https://doi.org/10.1016/00092509(65)850023
» https://doi.org/10.1016/00092509(65)850023  Haroon, S. S., Malik, T. N., Evaporation RateBased Water Cycle Algorithm for ShortTerm Hydrothermal Scheduling, Arabian Journal for Science and Engineering, 1, 116 (2016).
 Hu, X., Coello Coello, C.A., Huang, Z.A., New multiobjective evolutionary algorithm: neighborhood exploring evolution strategy, Engineering Optimization, 37, 351379 (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 multiobjective optimization algorithm based on differential evolution and neighborhood exploring evolution strategy. Journal of Artificial Intelligence and Soft Computing Research, 1, 112 (2011).
 Lobato, F.S., Steffen Jr., V., SilvaNeto, A. J., Solution of Singular Optimal Control Problems using the Improved Differential Evolution Algorithm. Journal of Artificial Intelligence and Soft Computing Research , 1, 2436 (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., Gradientbased Water Cycle Algorithm with evaporation rate applied to chaos suppression. Applied Soft Computing, 53, 420440 (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, 7579 (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, 470472 (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), 295303 (2008). https://doi.org/10.1590/S167858782008000400005
» https://doi.org/10.1590/S167858782008000400005  Ryoo, H. S., Sahinidis, B. P., Global Optimization of Nonconvex NLPs and MINLPs with Application in Process Design. Computers & Chemical Engineering, 19(5), 551566 (1995). https://doi.org/10.1016/00981354(94)000972
» https://doi.org/10.1016/00981354(94)000972  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, 5871 (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, 116 (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 multiobjective optimization problems, Soft Computing, 19(9), 25872603 (2015c). https://doi.org/10.1007/s0050001414244
» https://doi.org/10.1007/s0050001414244  Sadollah, A., Eskandar, H., Kim, J. H., Water cycle algorithm for solving constrained multiobjective problems, Applied Soft Computing , 27, 279298 (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, 103116 (2014).
 Vassiliadis, V., Computational Solution of Dynamic Optimization Problems with General DifferentialAlgebraic 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(34), 261275 (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, 125148 (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
JanMar 2019
History

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