On the Effectiveness of Nature-Inspired Metaheuristic Algorithms for Performing Phase Equilibrium Thermodynamic Calculations

The search for reliable and efficient global optimization algorithms for solving phase stability and phase equilibrium problems in applied thermodynamics is an ongoing area of research. In this study, we evaluated and compared the reliability and efficiency of eight selected nature-inspired metaheuristic algorithms for solving difficult phase stability and phase equilibrium problems. These algorithms are the cuckoo search (CS), intelligent firefly (IFA), bat (BA), artificial bee colony (ABC), MAKHA, a hybrid between monkey algorithm and krill herd algorithm, covariance matrix adaptation evolution strategy (CMAES), magnetic charged system search (MCSS), and bare bones particle swarm optimization (BBPSO). The results clearly showed that CS is the most reliable of all methods as it successfully solved all thermodynamic problems tested in this study. CS proved to be a promising nature-inspired optimization method to perform applied thermodynamic calculations for process design.


Introduction
Applied thermodynamic calculations in chemical engineering often involve the repeated solution of phase stability and phase equilibrium problems as their solutions are needed during the design of several equipment and separation processes. These problems can be formulated as minimization problems, for which the global minimum represents the required result. These calculations are challenging due to the high nonlinearity of thermodynamic models used to describe the equilibrium phases, the potential nonconvexity of the thermodynamic functions used as objective, and the presence of trivial solutions in the feasible search space. Thus, the solution of this type of problems via global optimization algorithms remains to be an active area of research. These problems generally feature local minima that are comparable to the global minimum, which accentuates the need for reliable global optimizers [1,2]. For example, the features of reactive phase equilibrium calculations increase the dimensionality and complexity of the optimization problem because the objective functions are required to satisfy the chemical equilibrium constraints [1,2].
The global stochastic optimization methods show high probabilities to locate the global minimum within reasonable computational costs, and thus they offer a desirable balance between reliability and efficiency for finding the global optimum solution. Moreover, stochastic methods do not require any assumptions for the optimization problem at hand, are more capable of addressing the nonlinearity and nonconvexity of the objective function, and are relatively easier to program and implement, among other advantages [3].
The application of stochastic global optimization methods for solving phase equilibrium thermodynamic problems has grown considerably during last years. To date, the most popular stochastic global optimization methods have been used and applied for solving phase equilibrium thermodynamic problems, for example, simulated annealing, genetic algorithms, tabu search, differential evolution, particle swarm optimization, and ant colony optimization (ACO) [4][5][6][7][8][9][10][11][12][13][14][15]. 2 The Scientific World Journal For example, a variant of ACO was tested in the global optimization of thermodynamic problems and was found to be robust in solving vapor-liquid equilibrium parameter estimation problems [4]. Zhu et al. [5] used an enhanced simulated annealing algorithm to solve multicomponent phase stability problems. Bonilla-Petriciolet and his coworkers compared different variants of PSO [6] and different variants of simulated annealing [14] for solving phase equilibrium problems. Repulsive particle swarm optimization was also studied by Rahman et al. [8]. Rangaiah and his co-workers studied the differential evolution [9,10], tabu search [11], and genetic algorithms [12] for solving phase stability and phase equilibrium problems.
The above studies have analyzed the capabilities and limitations of stochastic optimizers. But there exists no conclusive evaluation of those methods in comparison to one another for the solution of phase stability and phase equilibrium problems. Typically, each algorithm is introduced and compared with some of the other algorithms in a research publication. However, to the best of our knowledge, there exists no study that presents to the scientific community a ranking of the efficiency and reliability of those algorithms for the purpose of solving phase equilibrium and stability problems.
The aim of this study is provide a definitive ranking of the performance of a set of nature-inspired metaheuristic algorithms. To do so, we have selected eight of the most promising nature-inspired optimization methods based on the performance reported in the literature or obtained from our previous studies. These algorithms are cuckoo search (CS), intelligent firefly (IFA), bat (BA), artificial bee colony (ABC), monkey and krill herd hybrid (MAKHA), covariance matrix adaptation evolution strategy (CMAES), magnetic charged system search (MCSS), and bare bones particle swarm optimization (BBPSO). We systematically used those methods on some of the difficult phase stability and phase equilibrium problems reported in the literature and then analyzed their performance in terms of clear reliability and efficiency metrics.
The remainder of this paper is organized as follows. The eight optimization methods and the rationale for their selection are briefly presented in Section 2. A brief description of the phase stability and equilibrium problems is given in Section 3, including the implementation details of the eight algorithms. Section 4 presents the results and discussion of their performance in solving these thermodynamic calculations. Finally, the conclusions of this study are summarized in Section 5.

Selection and Description of the Nature-Inspired Metaheuristic Algorithms
Each of the eight selected metaheuristics is presented below. Only brief introductions are made here. Interested readers are referred to the primary sources of those algorithms for more information.
Cuckoo search (CS) is an optimization algorithm inspired by the obligate brood parasitism of some cuckoo species by laying their eggs in the nests of other host birds [16]. Intelligent firefly algorithm (IFA) [17] is a variant of firefly algorithm [18], a metaheuristic algorithm, inspired by the flashing behavior of fireflies to attract other fireflies. MAKHA is a hybrid between monkey algorithm (MA) [19], which is inspired by the simulation of the climbing processes of monkeys to find the highest mountaintop, and krill-herd algorithm (KHA) [20], which is based on the simulation of the herding behavior of krill individuals. Covariance matrix adaptation evolution strategy (CMAES) [21] is a stochastic and derivative free method for numerical optimization of nonlinear nonconvex problems. Artificial bee colony (ABC) [22] is an optimization algorithm based on the intelligent foraging behavior of honey bee swarm. Bat algorithm (BA) [23] is another bioinspired optimization algorithm based on the echolocation behavior of microbats with varying pulse rates of emission and loudness. Magnetic charged system search (MCSS) [24] is a variant of charged system search [25], which is based on the application of physics principles such as Coulomb law and Newtonian laws of mechanics to model how charged particles affect one another during their move towards the largest bodies. In MCSS, magnetic forces are also considered in addition to electrical forces. Finally, a variant of bare bones particle swarm optimization (BBPSO) [26] is based on the original particle swarm optimization [27], but without parameters and with the incorporation of mutation and crossover operators of DE to enhance the global search capability.
Since it was not possible to include all global stochastic optimization methods available in the literature for this comparative study, a screening process was performed to select the most promising ones. This process depended mainly on the results of solving phase stability and phase equilibrium problems using global optimization methods as reported in the literature. In several publications, limited comparisons were reported between some GSO methods. For example, CMAES was selected as it was shown to perform better than shuffled complex evolution in solving phase equilibrium and phase stability problems [28]; IFA performed better than FA in general [17], CS better than integrated differential evolution [29], MCSS better than CSS for phase equilibrium and phase stability problems [30], and BBPSO better than PSO [26]. In addition, our preliminary calculations showed that MAKHA performed better than MA and KHA, and ABC and BA performed better than FA.
One approach to solving phase stability and phase equilibrium problems is to start the optimization process with a stochastic global optimizer, as the methods studied in this work. Once a certain stopping criterion is satisfied, we follow with a local optimizer, such as sequential quadratic programming, to close down to the minimum within the vicinity of the best value found by the global optimizer. This approach has been proven successful in previous studies [28][29][30] and it would complement any of the methods studied above. However, we restricted this study to the performance of the stochastic global optimizers without the use of a local optimizer to focus on the strength and weakness of the studied methods free from any artificial enhancement of their results.
The Scientific World Journal 3

Objective Functions.
In this study, the phase stability and equilibrium problems are stated as a global optimization problem. Therefore, the global optimization problem to be solved is as follows: minimize (X) with respect to decision variables: X = ( 1 , . . . , ). The upper and lower bounds of these variables are ( 1 max , . . . , max ) and ( 1 min , . . . , min ), respectively.
The phase stability, phase equilibrium, and reactive phase equilibrium calculations for testing the performance of global optimization methods are explained briefly in Table 1, which shows the problem formulation, objective function, decision variables, and constraints used for those thermodynamic calculations. Specifically, the phase stability analysis was performed using the global minimization of the tangent plane distance function (TPDF) [31], while the global optimization of the Gibbs free energy was used for phase equilibrium calculations with or without chemical reactions [2]. The mathematical formulation for phase stability and phase equilibrium calculations for nonreactive systems is an unconstrained minimization of the objective function, while the constrained Gibbs free energy minimization in reactive systems was performed using the penalty function method according to the approach reported by Bonilla-Petriciolet et al. [1]. For interested readers, several references provide a detailed description of these thermodynamic calculations [1,2,4,10,12].
Previous work reported the evaluation of global optimization methods for solving twenty-four problems [4,28,30]. In this work, we focused on the nine most difficult ones. The basis for the selection was the relatively lower success rates that optimization methods obtained when solving them in the previous studies. These problems are presented in Table 2.

Details of Numerical Implementation and Performance
Metrics Used for Testing the Algorithms. All thermodynamic problems and the different optimization algorithms were coded in the MATLAB technical computing environment. The codes for CS and BA were obtained from MATLAB file exchange server as uploaded by their developers and used without change. The code for IFA was developed by the authors through minor modifications of the FA code that was obtained from the MATLAB file exchange server as well. The codes for CMAES and ABC were obtained from the developers' web sites and used without change. The code for MCSS was written by the authors based on the developer's published work [24,25]. MAKHA was developed and coded by the authors. The code for BBPSO was obtained from its developer [26]. Each problem was solved 30 times independently and with different random initial seeds to determine the reliability of the optimization algorithms. Calculations were performed for a certain number of iterations and then stopped. This maximum value for the number of iterations was different for different algorithms. The maximum values were selected to give the same number of function evaluations at the end of the run. Table 3 shows the values selected for the parameters of the eight optimization algorithms, which were determined using preliminary calculations.
The methods were evaluated according to the reliability and efficiency for finding the global optimum. The efficiency is determined by recording the number of function evaluations NFE for each optimization algorithm, where a low value of NFE means a higher efficiency. Note that NFE is an unbiased indicator of the computational costs required by a certain algorithm and is independent of the host hardware. In previous studies [1,4,6,26,28,30], reliability was measured by the success rate at certain number of iterations. The success rate is defined as the ratio of number of runs in which the global minimum was attained within a tolerance at this iteration number to the total number of runs. In this work, we present a different reliability metric: a plot of the average best value against the number of function evaluations. The best values are averaged over all the runs and plotted against NFE, which is calculated at each iteration. Since the NFE needed for each iteration differs amongst the optimization methods, the plot of average best value against NFE is a better indication of reliability versus efficiency of the optimization method.
For a comparative evaluation of the global optimization methods, we have employed performance profile (PP) reported by Dolan and Moré [32], who introduced PP as a tool for evaluating and comparing the performance of optimization software. In particular, PP has been proposed to represent compactly and comprehensively the data collected from a set of solvers for a specified performance metric such as the computing time or the number of function evaluations. The PP plot allows visualization of the expected performance differences among several solvers and comparing the quality of their solutions by eliminating the bias of failures obtained in a small number of problems.
Consider solvers (i.e., optimization methods) to be tested over a set of problems. For each problem and solver , the performance metric must be defined. In our study, reliability of the stochastic method in accurately finding the global minimum of the objective function is considered as the principal goal, and hence the reliability performance metric is defined as where * is the known global optimum of the objective function and calc is the mean value of that objective function calculated by the metaheuristic over several runs. We have also used another performance metric for the evaluation of the efficiency of the method in obtaining the global minimum. This metric is the minimum number of NFE needed to reach with 10 −5 of the global minimum.
For the performance metric of interest, the performance ratio, , is used to compare the performance on problem by solver with the best performance by any solver on this problem. This performance ratio is given by 4 The Scientific World Journal    The Scientific World Journal The value of is 1 for the solver that performs the best on a specific problem . To obtain an overall assessment of the performance of solvers on problems, the following cumulative function for is used: where ( ) is the fraction of the total number of problems, for which solver has a performance ratio within a factor of of the best possible ratio. The PP of a solver is a plot of ( ) versus ; it is a nondecreasing, piecewise constant function, continuous from the right at each of the breakpoints [32]. To identify the best solver, it is only necessary to compare the values of ( ) for all solvers and to select the highest one, which is the probability that a specific solver will "win" over the rest of solvers used.
In our case, one PP plot compares how accurately the stochastic methods can find the global optimum value relative  to one another, and so the term "win" refers to the stochastic method that provides the most accurate value of the global minimum in the benchmark problems used. The other PP plot compares how fast the stochastic methods can find the global minimum with a tolerance level of 10 −5 , so the term "win", in this case, refers to the method that reaches the solution fastest for the problems used.

Results and Discussion
The results are presented in three different ways. For each problem, the mean best values are plotted versus NFE for each of the eight algorithms. These plots are found in Figures 1-9. The minimum NFE required to reach a certain tolerance from the known global minimum for each problem was calculated and presented in Table 4. The performance profiles for the reliability and efficiency metrics are shown in Figures 10 and  11, respectively. A detailed discussion of the results follows.
The Scientific World Journal

Phase Stability Problems.
Problem T7 is a nine-variable phase-stability problem that is extremely difficult to solve. The means of the minimum values obtained by all methods were not close enough to the global minimum except for CS. As shown in Figure 1 and Table 4, ABC and MCSS were able to get to within 10 −3 of the global minimum. On the other hand, CS was able to find the global minimum down to a tolerance of 10 −7 . To reach the global minimum within a tolerance of 10 −5 , it required 109280 function evaluations. Problem T8 is also a difficult phase stability problem. Figure 2 shows how all problems were able to reach values close to the global optimum. However, close analysis at the vicinity of the global minimum, as depicted in the inset of Figure 2, at the level of 10 −5 , revealed that MAKHA and BA failed to find the global minimum up to the end of the runs. CMAES was the most efficient as it converged to the global minimum in the least NFE by at least one order of magnitude. None of the methods was able to reach within 10 −6 of the global minimum, as shown in Table 4.
Problem T9 is the last of the three phase stability problems. Even though, MAKHA was quite fast in approaching the global minimum, as depicted in Figure 3, it failed at converging to within 10 −5 of the global minimum. IFA was also not able to find the global minimum. CMAES was the most efficient method in getting down to 10 −5 distance from the global minimum but was not able to get any closer. CS, again, was the only method to converge reliably down to 10 −7 of the global minimum.
For the phase stability problems, CS is clearly the most reliable method. It may not be as efficient in its initial approach to the global minimum as other methods such as BA or CMAES, but it outperforms the rest in terms of finding the global minimum. An open area of development for CS would be to make it more efficient via hybridization with some of the other methods in their initial approach to the global minimum.

Metaheuristic
Phase equilibrium thermodynamic problem  T7  T8  T9  G4  G6  G7  G8  R4  R7 4.2. Phase Equilibrium Problems. Problem G4 is a twovariable phase equilibrium problem that is relatively easy to solve. However, CMAES seemed to have been trapped in a local minimum and was unable to find its global minimum, within a tolerance of 10 −5 , as shown in Figure 4. IFA did slightly better than CMAES, but was unable to reach the global minimum within a tolerance of 10 −6 . MAKHA was the most efficient in finding the global minimum within 10 −6 and 10 −7 , with BBPSO and CS performing quite well.
Despite the fact that CMAES was not able to solve problem G4, it was superior in solving problem G6. With only 101 NFE, CMAES reached down to 10 −6 of the global minimum, as is shown in Figure 5. All methods converged to 10 −6 from the global minimum, but only CMAES, CS, and MCSS converged to 10 −7 , with CMAES being ten times more efficient. This convergence pattern was repeated in problem G7. Only CMAES and CS solved the problem down to the 10 −6 and 10 −7 levels, with CMAES being one order of magnitude 10 The Scientific World Journal more efficient, as is clear in Figure 6 and Table 4. MAKHA, BA, and BBPSO were not able to converge at the 10 −5 level.
Problem G8 was successfully solved at the 10 −5 level by IFA, CMAES, ABC, BA, CS, and BBPSO, as shown in Figure 7. Only CMAES and CS solved the problem down to the 10 −7 levels, with CMAES being one order of magnitude more efficient. In fact, CMAES was quite efficient at all tolerance levels, as shown by the NFE numbers in Table 4.
The convergence profiles of the four phase equilibrium problems (G4, G6, G7, and G8) indicated that CS is the most reliable of all algorithms as it was the only one to be able to solve all problems down to the 10 −7 tolerance level. CMAES was the most efficient as it required one order of magnitude less NFE to solve three of the four problems down to the same tolerance level. However, CMAES failed to solve the two-variable problem that was successfully solved by all other methods, except IFA, down to the 10 −7 level.

Reactive Phase Equilibrium Problems.
Regardless of the number of variables, the reactive phase equilibrium problems are more difficult than the nonreactive phase equilibrium problems because the chemical reaction equilibria constraints must be satisfied. Problem R4, see Figure 8, was successfully solved down to the 10 −5 tolerance level by CS, which was also able to converge to the global minimum at the 10 −6 and 10 −7 levels. MAKHA, CMAES, BA, MCSS, and BBPSO were not able to arrive even at a level of 10 −3 from the global minimum. Similarly, CMAES and BA were not able to reach the 10 −3 level for Problem R7. However, MAKHA, IFA, CS, and BBPSO converged down to 10 −7 distance from the global minimum, with IFA being the most efficient down to the 10 −5 level and BBPSO at the 10 −6 and 10 −7 levels.
The complete failure of CMAES to solve reactive phase equilibrium problems is remarkable. CMAES functions extremely well in certain types of problems and extremely bad in others. On the other hand, CS solved the reactive phase equilibrium problems just as it reliably solved all other problems in this study. Since CS uses Lévy walk, instead of random walk, in its global search, it can explore the search space more efficiently and avoid entrapment in local minima, as was demonstrated by our results. However, CS requires significantly large NFE to allow it to converge to the global minimum. Any attempt to improve CS performance should target its slow convergence behavior.
Our results are summarized in the PP plots of Figures  9 and 10. The reliability ranking, as extracted from the reliability PP plot of Figure 9, is as follows. CS is the most reliable, followed by CMAES, BBPSO, and MCSS, on the second level. The third level contains MAKHA, ABC, IFA, and BA, in that order. The efficiency ranking starts with CMAES, BBPSO, and ABC. The second level contains CS and IFA. The third level contains BA, MAKHA, and MCSS.

Conclusions
In this study, we have selected eight promising natureinspired metaheuristic algorithms for the solution of nine difficult phase stability and phase equilibrium problems. These thermodynamic problems were systematically solved by the different metaheuristics and the results were tracked and compared. The results clearly show that CS is the most reliable of all tested optimization methods as it successfully solved all problems down to the 10 −5 tolerance from the global minima. Any attempt to improve the performance of CS should target its slow convergence behavior. Recently developed CS variants [33] could provide more efficient performance for the solution of phase stability and phase equilibrium problems. These variants could be evaluated in a future study in an attempt to find the most reliable and efficient algorithm for this application. On the other hand, CMAES was the most efficient in finding the solution for the problems it was able to solve. However, it was not able to converge to the global minimum for some of the tested thermodynamic problems.

Nomenclature
: Parameter used in BA ABC: Artificial bee colony : Parameter used in MAKHA algorithm BA: Bat algorithm BBPSO: Bare bones particle swarm optimization : Parameter used in MAKHA algorithm CMAES: Covariance matrix adaptation evolution strategy CMCR: Parameter used in MCSS CMS: Parameter used in MCSS CS: Cuckoo search : Empirical and experimental constants-used in MAKHA : Dimension of the problem, Parameter used in MAKHA Algorithm max : Maximum diffusion speed-used in MAKHA , obj : Objective function FA: Firefly algorithm : Gibbs free energy g * : Global minimum : Ac o u n t e r , : Index of the component, or index used in an algorithm IFA: Intelligent firefly algorithm eq : Equilibriumconstant : Dimension of the problem, that is, number of variables MAKHA: Monkey Algorithm-Krill Herd Algorithm hybrid MCSS: Magnetic charged system search : The Numberofsolvers NV: Dimension of the problem, that is, number of variables : Parameter used in CS PP: Performance profile : The half-range of boundaries between the lower boundary and the upper boundary of the decision variables ( )-used in MAKHA : Parameter used in BA : Theperformanceratio TPDF: Tangent plane distance function : Performancemetric : Foraging speed-used in MAKHA : Inertiaweight-usedinMAKHA : Decisionvariables : Decisionvariables : Decisionvariables : Decisionvariables : Numberofcomponents.

Greek Letters
: The simulating value of : Thecounterof points max : The maximum assumed value of : Transformed decision variables used instead of mole fractions : Vector of random numbers in FA : Fugacity coefficient, Fraction of top fireflies to be utilized in the move-Parameter used in IFA : Activity coefficient, Parameter used in IFA : Numberofphases : Chemicalpotential : Parameter used in IFA : Parameter used in IFA : The cumulative probabilistic function of and the fraction of the total number of problems : Parameter used in CMAES.

Subscripts
: Feed : Index for the components in the mixture min: Minimum value : Initialvalueoftheparameter : At composition : At composition .