1. Introduction {#sec1}
===============

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 \[[@B1], [@B2]\]. 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 \[[@B1], [@B2]\].

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 \[[@B3]\].

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) \[[@B4]--[@B15]\]. 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 \[[@B4]\]. Zhu et al. \[[@B5]\] used an enhanced simulated annealing algorithm to solve multicomponent phase stability problems. Bonilla-Petriciolet and his coworkers compared different variants of PSO \[[@B6]\] and different variants of simulated annealing \[[@B14]\] for solving phase equilibrium problems. Repulsive particle swarm optimization was also studied by Rahman et al. \[[@B8]\]. Rangaiah and his co-workers studied the differential evolution \[[@B9], [@B10]\], tabu search \[[@B11]\], and genetic algorithms \[[@B12]\] 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](#sec2){ref-type="sec"}. A brief description of the phase stability and equilibrium problems is given in [Section 3](#sec3){ref-type="sec"}, including the implementation details of the eight algorithms. [Section 4](#sec4){ref-type="sec"} presents the results and discussion of their performance in solving these thermodynamic calculations. Finally, the conclusions of this study are summarized in [Section 5](#sec5){ref-type="sec"}.

2. Selection and Description of the Nature-Inspired Metaheuristic Algorithms {#sec2}
============================================================================

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 \[[@B16]\]. Intelligent firefly algorithm (IFA) \[[@B17]\] is a variant of firefly algorithm \[[@B18]\], a metaheuristic algorithm, inspired by the flashing behavior of fireflies to attract other fireflies. MAKHA is a hybrid between monkey algorithm (MA) \[[@B19]\], which is inspired by the simulation of the climbing processes of monkeys to find the highest mountaintop, and krill-herd algorithm (KHA) \[[@B20]\], which is based on the simulation of the herding behavior of krill individuals. Covariance matrix adaptation evolution strategy (CMAES) \[[@B21]\] is a stochastic and derivative free method for numerical optimization of nonlinear nonconvex problems. Artificial bee colony (ABC) \[[@B22]\] is an optimization algorithm based on the intelligent foraging behavior of honey bee swarm. Bat algorithm (BA) \[[@B23]\] 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) \[[@B24]\] is a variant of charged system search \[[@B25]\], 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) \[[@B26]\] is based on the original particle swarm optimization \[[@B27]\], 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 \[[@B28]\]; IFA performed better than FA in general \[[@B17]\], CS better than integrated differential evolution \[[@B29]\], MCSS better than CSS for phase equilibrium and phase stability problems \[[@B30]\], and BBPSO better than PSO \[[@B26]\]. 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 \[[@B28]--[@B30]\] 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.

3. Description of Phase Stability and Phase Equilibrium Problems Used for the Evaluation {#sec3}
========================================================================================

3.1. Objective Functions {#sec3.1}
------------------------

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 *F*(**X**) with respect to *D* decision variables: **X** = (*X* ^1^, ..., *X* ^*D*^). The upper and lower bounds of these variables are (*X* ~max⁡~ ^1^,..., *X* ~max⁡~ ^*D*^) and (*X* ~min⁡~ ^1^,..., *X* ~min⁡~ ^*D*^), 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](#tab1){ref-type="table"}, 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) \[[@B31]\], while the global optimization of the Gibbs free energy was used for phase equilibrium calculations with or without chemical reactions \[[@B2]\]. 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. \[[@B1]\]. For interested readers, several references provide a detailed description of these thermodynamic calculations \[[@B1], [@B2], [@B4], [@B10], [@B12]\].

Previous work reported the evaluation of global optimization methods for solving twenty-four problems \[[@B4], [@B28], [@B30]\]. 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](#tab2){ref-type="table"}.

3.2. Details of Numerical Implementation and Performance Metrics Used for Testing the Algorithms {#sec3.2}
------------------------------------------------------------------------------------------------

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 \[[@B24], [@B25]\]. MAKHA was developed and coded by the authors. The code for BBPSO was obtained from its developer \[[@B26]\]. 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](#tab3){ref-type="table"} 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 \[[@B1], [@B4], [@B6], [@B26], [@B28], [@B30]\], 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é \[[@B32]\], 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 *n* ~*s*~ solvers (i.e., optimization methods) to be tested over a set of *n* ~*p*~ problems. For each problem *p* and solver *s*, the performance metric *t* ~*ps*~ 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 $$\begin{matrix}
{t_{ps} = f_{\text{calc}} - f^{\ast},} \\
\end{matrix}$$ where *f*\* is the known global optimum of the objective function and *f* ~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, *r* ~*ps*~, is used to compare the performance on problem *p* by solver *s* with the best performance by any solver on this problem. This performance ratio is given by $$\begin{matrix}
{r_{ps} = \frac{t_{ps}}{{\min}\left\{ {t_{ps}:1 \leq s \leq n_{s}} \right\}}.} \\
\end{matrix}$$

The value of *r* ~*ps*~ is 1 for the solver that performs the best on a specific problem *p*. To obtain an overall assessment of the performance of solvers on *n* ~*p*~ problems, the following cumulative function for *r* ~*ps*~ is used: $$\begin{matrix}
{\rho_{s}\left( \zeta \right) = \frac{1}{n_{p}}\text{size}\left\{ {p:r_{ps} \leq \zeta} \right\},} \\
\end{matrix}$$ where *ρ* ~*s*~(*ζ*) is the fraction of the total number of problems, for which solver *s* has a performance ratio *r* ~*ps*~ within a factor of *ζ* of the best possible ratio. The PP of a solver is a plot of *ρ* ~*s*~(*ζ*) versus *ζ*; it is a nondecreasing, piecewise constant function, continuous from the right at each of the breakpoints \[[@B32]\]. To identify the best solver, it is only necessary to compare the values of *ρ* ~*s*~(*ζ*) 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.

4. Results and Discussion {#sec4}
=========================

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](#fig1){ref-type="fig"}--[9](#fig9){ref-type="fig"}. The minimum NFE required to reach a certain tolerance from the known global minimum for each problem was calculated and presented in [Table 4](#tab4){ref-type="table"}. The performance profiles for the reliability and efficiency metrics are shown in Figures [10](#fig10){ref-type="fig"} and [11](#fig11){ref-type="fig"}, respectively. A detailed discussion of the results follows.

4.1. Phase Stability Problems {#sec4.1}
-----------------------------

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](#fig1){ref-type="fig"} and [Table 4](#tab4){ref-type="table"}, 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](#fig2){ref-type="fig"} 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](#fig2){ref-type="fig"}, 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](#tab4){ref-type="table"}.

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](#fig3){ref-type="fig"}, 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.

4.2. Phase Equilibrium Problems {#sec4.2}
-------------------------------

Problem G4 is a two-variable 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](#fig4){ref-type="fig"}. 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](#fig5){ref-type="fig"}. 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 more efficient, as is clear in [Figure 6](#fig6){ref-type="fig"} and [Table 4](#tab4){ref-type="table"}. 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](#fig7){ref-type="fig"}. 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](#tab4){ref-type="table"}.

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.

4.3. Reactive Phase Equilibrium Problems {#sec4.3}
----------------------------------------

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](#fig8){ref-type="fig"}, 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](#fig9){ref-type="fig"} and [10](#fig10){ref-type="fig"}. The reliability ranking, as extracted from the reliability PP plot of [Figure 9](#fig9){ref-type="fig"}, 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.

5. Conclusions {#sec5}
==============

In this study, we have selected eight promising nature-inspired 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 \[[@B33]\] 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.

*A*:

:   Parameter used in BA

ABC:

:   Artificial bee colony

*B*:

:   Parameter used in MAKHA algorithm

BA:

:   Bat algorithm

BBPSO:

:   Bare bones particle swarm optimization

*C*:

:   Parameter used in MAKHA algorithm

CMAES:

:   Covariance matrix adaptation evolution strategy

CMCR:

:   Parameter used in MCSS

CMS:

:   Parameter used in MCSS

CS:

:   Cuckoo search

*C*~*t*~:

:   Empirical and experimental constants---used in MAKHA

*D*:

:   Dimension of the problem, Parameter used in MAKHA Algorithm

*D*~max⁡~:

:   Maximum diffusion speed---used in MAKHA

*F*, *F*~obj~:

:   Objective function

FA:

:   Firefly algorithm

*G*:

:   Gibbs free energy

**g**\*:

:   Global minimum

*I*:

:   A counter

*i*, *j*:

:   Index of the component, or index used in an algorithm

IFA:

:   Intelligent firefly algorithm

*K*~eq~:

:   Equilibrium constant

*m*:

:   Dimension of the problem, that is, number of variables

MAKHA:

:   Monkey Algorithm-Krill Herd Algorithm hybrid

MCSS:

:   Magnetic charged system search

*N*:

:   The number of iterations (criterion maximum number)

*n*:

:   Population number

NFE:

:   Number of function evaluations

*N*~max⁡~:

:   Maximum induced speed

*N*~*P*~:

:   Population size (number of points)

*n*~*p*~:

:   Number of problems

**n**~ref~:

:   Column vector of moles of each of the reference components

*n*~*s*~:

:   Number of solvers

NV:

:   Dimension of the problem, that is, number of variables

*p*:

:   Parameter used in CS

PP:

:   Performance profile

*R*:

:   The half-range of boundaries between the lower boundary and the upper boundary of the decision variables (*X*)---used in MAKHA

*r*:

:   Parameter used in BA

*r*~*ps*~:

:   The performance ratio

TPDF:

:   Tangent plane distance function

*t*~*ps*~:

:   Performance metric

*V*~*f*~:

:   Foraging speed---used in MAKHA

*w*~*f*~:

:   Inertia weight---used in MAKHA

*X*:

:   Decision variables

*x*:

:   Decision variables

*Y*:

:   Decision variables

*y*:

:   Decision variables

*c*:

:   Number of components.

```{=html}
<!-- -->
```

*ζ*:

:   The simulating value of *r* ~*ps*~

*ς*:

:   The counter of *ρ* points

*ζ*~max⁡~:

:   The maximum assumed value of *r* ~*ps*~

*β*:

:   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

*π*:

:   Number of phases

*μ*:

:   Chemical potential

*α*~*o*~:

:   Parameter used in IFA

*β*~*o*~:

:   Parameter used in IFA

*ρ*:

:   The cumulative probabilistic function of *r* ~*ps*~ and the fraction of the total number of problems

*σ*:

:   Parameter used in CMAES.

```{=html}
<!-- -->
```

*F*:

:   Feed

*I*:

:   Index for the components in the mixture

min:

:   Minimum value

*O*:

:   Initial value of the parameter

*y*:

:   At composition *y*

*z*:

:   At composition *z*.

```{=html}
<!-- -->
```

0:

:   Pure component.

Conflict of Interests
=====================

The authors of this paper do not have a direct financial relation that might lead to a conflict of interests with regard to this paper.

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem T7.](TSWJ2014-374510.001){#fig1}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem T8.](TSWJ2014-374510.002){#fig2}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem T9.](TSWJ2014-374510.003){#fig3}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem G4.](TSWJ2014-374510.004){#fig4}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem G6.](TSWJ2014-374510.005){#fig5}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem G7.](TSWJ2014-374510.006){#fig6}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem G8.](TSWJ2014-374510.007){#fig7}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem RG4.](TSWJ2014-374510.008){#fig8}

![The evolution of the mean best value calculated via the eight metaheuristics versus NFE for problem RG7.](TSWJ2014-374510.009){#fig9}

![Performance profile (PP) of the reliability metric of the eightmetaheuristics for the 9 phase stability and equilibrium problems.](TSWJ2014-374510.010){#fig10}

![Performance profile (PP) of the efficiency metric of the eight metaheuristics for the 9 phase stability and equilibrium problems.](TSWJ2014-374510.011){#fig11}

###### 

Description of thermodynamic functions and optimization problems for phase stability analysis and equilibrium calculations in reactive and nonreactive systems.

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Calculation                               Description                                                                                                                                                                                                   Thermodynamic function                                                                                                                                                                                                                                                                                                                                                                         Optimization problem
  ----------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Phase stability analysis                  It involves the determination of whether a system will remain in one phase at the given conditions or split into two or more phases                                                                           Tangent plane distance function\                                                                                                                                                                                                                                                                                                                                                               $\begin{matrix}
                                                                                                                                                                                                                                                          $\text{TPDF} = {\sum\limits_{i = 1}^{c}{y_{i}\left( {\left. \mu_{i} \right|_{y} - \left. \mu_{i} \right|_{z}} \right)}}$  \                                                                                                                                                                                                                                                                    {\min\limits_{\beta}\text{TPDF\quad\quad\quad\quad\quad}} \\
                                                                                                                                                                                                                                                          where *c* is the number of components of the mixture and μ~*i*~\|~*y*~ and μ~*i*~\|~*z*~ are the chemical potentials calculated at trial composition *y* and feed composition *z*                                                                                                                                                                                                              {0 \leq \beta_{i} \leq 1\quad i = 1,\ldots,c.} \\
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         \end{matrix}$ \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         The decision variables are *β* ~*i*~ ∈ (0,1) using the following relationships: \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         *n* ~*iy*~ = β~*i*~ *z* ~*i*~ *n* ~*F*~  *i* = 1,..., *c*  \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         $y_{i} = \frac{n_{iy}}{\sum_{j = 1}^{c}n_{iy}}\quad i = 1,\ldots,c$  \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         where *n* ~*iy*~ are the mole numbers of component *i* in phase *y* and *n* ~*F*~ is the total moles in the mixture under analysis

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         

  Phase equilibrium calculation             It involves the determination of the number, type, and composition of the phases at equilibrium at the given operating conditions                                                                             Gibbs free energy of mixing (*g*) \                                                                                                                                                                                                                                                                                                                                                            $\begin{matrix}
                                                                                                                                                                                                                                                          $\begin{matrix}                                                                                                                                                                                                                                                                                                                                                                                {\min\limits_{\beta}{g\text{\quad\quad\quad\quad\quad}}} \\
                                                                                                                                                                                                                                                          {g = \sum\limits_{j = 1}^{\pi}\sum\limits_{i = 1}^{c}{n_{ij}{\ln{(x_{ij}\gamma_{ij})}}}} \\                                                                                                                                                                                                                                                                                                    {0 \leq \beta_{ij} \leq 1\text{\quad\quad\quad}} \\
                                                                                                                                                                                                                                                          {\quad   = {\sum\limits_{j = 1}^{\pi}{\sum\limits_{i = 1}^{c}{n_{ij}{\ln\left( \frac{x_{ij}{\hat{\varphi}}_{ij}}{\varphi_{i}} \right)}}}}} \\                                                                                                                                                                                                                                                  {i = 1,\ldots,c\text{\quad\quad\quad}} \\
                                                                                                                                                                                                                                                          \end{matrix}$   \                                                                                                                                                                                                                                                                                                                                                                              {j = 1,\ldots,\pi - 1.\text{\quad\quad}} \\
                                                                                                                                                                                                                                                          where *π* is the number of phases at equilibrium and *θ* ~*ij*~ denotes the composition (i.e., *x* or *n*) or thermodynamic property of component *i* in phase *j*                                                                                                                                                                                                                             \end{matrix}$ \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         The decision variables are *β* ~*ij*~ ∈ (0,1) using the following relationships:\
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         *n* ~*i*1~ = β~*i*1~ *z* ~*i*~ *n* ~*F*~  *i* = 1, \..., *c* \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         $\begin{matrix}
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         {n_{ij} = \beta_{ij}\left( {z_{i}n_{F} - \sum\limits_{m = 1}^{j - 1}}n_{im} \right)\text{\quad  }} \\
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         {i = 1,\ldots,c;\quad j = 2,\ldots,\pi - 1} \\
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         \end{matrix}$ \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         $n_{i\pi} = z_{i}n_{F} - \sum\limits_{m = 1}^{\pi - 1}{n_{im}\quad}i = 1,...,c$

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         

  Reactive phase equilibrium Calculations   It involves the determination of the number, type and composition of the phases at equilibrium at the given operating conditions and subject to element/mass balances and chemical equilibrium constraints.   Gibbs free energy of mixing defined using reaction equilibrium constants \[[@B2]\] \                                                                                                                                                                                                                                                                                                           $\min\limits_{n_{ij}}G_{K}$ \
                                                                                                                                                                                                                                                          $G_{K} = g - \sum\limits_{j = 1}^{\pi}{\ln{\mathbf{K}_{\text{eq}}\mathbf{N}^{- 1}\mathbf{n}_{\text{ref},j}}}$ \                                                                                                                                                                                                                                                                                subject to \
                                                                                                                                                                                                                                                          where *g* is the Gibbs free energy of mixing, ln⁡**K** ~eq~ is a row vector of logarithms of chemical equilibrium constants for *r* independent reactions, **N** is an invertible, square matrix formed from the stoichiometric coefficients of a set of reference components chosen from the *r* reactions, and **n** ~ref~ is a column vector of moles of each of the reference components   $\sum\limits_{j = 1}^{\pi}\left( {n_{ij} - \mathbf{v}_{i}\mathbf{N}^{- 1}\mathbf{n}_{\text{ref},j}} \right) = n_{iF} - \mathbf{v}_{i}\mathbf{N}^{- 1}\mathbf{n}_{\text{ref},F}$ \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         *i* = 1,..., *c* − *r*  \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         *n* ~*ij*~ \> 0 *i* = 1,..., c *j* = 1,..., π  \
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         where *n* ~*i*,*F*~ is the initial moles of component *i* in the feed, **v** ~*i*~ is the row vector (of dimension *r*) of stoichiometric coefficients of component *i* in *r* reactions, and *n* ~*ij*~ is the number of moles of component *i* in phase *j*. The constrained global optimization problem can be solved by minimizing *G* ~*K*~ with respect to *c*(*π* − 1) + *r* decision variables *n* ~*ij*~. In this formulation, the mass balance equations are rearranged to reduce the number of decision variables of the optimization problem and to eliminate equality constraints
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

###### 

Details of the phase stability, phase equilibrium, and reactive phase equilibrium problems used in this study.

  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Code   System                                                                   Feed conditions                                                                                              Thermodynamic models                                                     Global optimum
  ------ ------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------ ----------------
  T7     C~1~ + C~2~ + C~3~ + C~4~ + C~5~ + C~6~ + C~7--16~ + C~17+~              *n* ~*F*~ = (0.7212, 0.09205, 0.04455, 0.03123, 0.01273, 0.01361, 0.07215, 0.01248) at 353 K and 38500 kPa   Phase stability problem with SRK EoS with classical mixing rules         −0.002688

                                                                                                                                                                                                                                                                        

  T8     C~1~ + C~2~ + C~3~ + *i*C~4~ + C~4~ + *i*C~5~ + C~5~ + C~6~ + *i*C~15~   *n* ~*F*~ = (0.614, 0.10259, 0.04985, 0.008989, 0.02116, 0.00722, 0.01187, 0.01435, 0.16998)\                Phase stability problem with SRK EoS with classical mixing rules         −1.486205
                                                                                  at 314 K and 2010.288 kPa                                                                                                                                                             

                                                                                                                                                                                                                                                                        

  T9     C~1~ + C~2~ + C~3~ + C~4~ + C~5~ + C~6~ + C~7~ + C~8~ + C~9~ + C~10~     *n* ~*F*~ = (0.6436, 0.0752, 0.0474, 0.0412, 0.0297, 0.0138, 0.0303, 0.0371, 0.0415, 0.0402) \               Phase stability problem with SRK EoS with classical mixing rules         −0.0000205
                                                                                  at 435.35 K and 19150 kPa                                                                                                                                                             

                                                                                                                                                                                                                                                                        

  G4     C~1~ + H~2~S                                                             *n* ~*F*~ = (0.9813, 0.0187)  \                                                                              Phase equilibrium problem with SRK EoS with classical mixing rules       −0.019892
                                                                                  at 190 K and 4053 kPa                                                                                                                                                                 

                                                                                                                                                                                                                                                                        

  G6     C~2~ + C~3~ + C~4~ + C~5~ + C~6~                                         *n* ~*F*~ = (0.401, 0.293, 0.199, 0.0707, 0.0363)  \                                                         Phase equilibrium problem with SRK EoS with classical mixing rules       −1.183653
                                                                                  at 390 K and 5583 kPa                                                                                                                                                                 

                                                                                                                                                                                                                                                                        

  G7     C~1~ + C~2~ + C~3~ + C~4~ + C~5~ + C~6~ + C~7--16~ + C~17+~              *n* ~*F*~ = (0.7212, 0.09205, 0.04455, 0.03123, 0.01273, 0.01361, 0.07215, 0.01248) at 353 K and 38500 kPa   Phase equilibrium problem with SRK EoS with classical mixing rules       −0.838783

                                                                                                                                                                                                                                                                        

  G8     C~1~ + C~2~ + C~3~ + *i*C~4~ + C~4~ + *i*C~5~ + C~5~ + C~6~ + *i*C~15~   *n* ~*F*~ = (0.614, 0.10259, 0.04985, 0.008989, 0.02116, 0.00722, 0.01187, 0.01435, 0.16998)\                Phase equilibrium problem with SRK EoS with classical mixing rules       −0.769772
                                                                                  at 314 K and 2010.288 kPa                                                                                                                                                             

                                                                                                                                                                                                                                                                        

  R4     A1 + A2 *↔* A3 + A4 \                                                    *n* ~*F*~ = (0.3, 0.4, 0.3, 0.0) at 298.15 K and 101.325 kPa                                                 Reactive phase equilibrium problem with UNIQUAC model and ideal gas: \   −1.10630
         (1) Acetic acid\                                                                                                                                                                      ln⁡*K* ~eq~ = $\frac{450}{T}$ + 0.8                                      
         (2) n-Butanol\                                                                                                                                                                                                                                                 
         (3) Water\                                                                                                                                                                                                                                                     
         (4) n-Butyl acetate                                                                                                                                                                                                                                            

                                                                                                                                                                                                                                                                        

  R7     A1 + A2 *↔* A3                                                           *n* ~*F*~ = (0.52, 0.48, 0.0) at 323.15 K and 101.325 kPa                                                    Reactive phase equilibrium problem with Margules solution model:  \      −0.653756
                                                                                                                                                                                               *K* ~eq~ = 3.5                                                           
  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

###### 

Selected values of the parameters used in the implementation of the eight nature-inspired metaheuristic algorithms.

  Metaheuristic                    Parameter            Selected value
  -------------------------------- -------------------- ----------------
  MAKHA                            *n*                  40*D* 
  *B*                              0.5                  
  *C*                              −0.1                 
  *D*                              0.1                  
  *D* ~max⁡~                       0                    
  *C* ~*t*~                        0.5                  
  *V* ~*f*~                        0.2                  
  *W* ~*f*~                        0.1                  
                                                        
  IFA                              *α* ~*o*~            0.5 
  *β* ~min⁡~                       0.2                  
  *γ*                              1                    
  *n*                              20*D*                
  *ϕ*                              0.05                 
                                                        
  CMAES                            *σ*                  0.2
  *n*                              20*D*                
                                                        
  ABC                              *n*                  20*D* 
  Food number                      *n*/2                
  limit                            100                  
                                                        
  BA                               *n*                  20*D* 
  *A*                              0.25                 
  *r*                              0.5                  
                                                        
  MCSS                             CMCR                 0.95 
  PAR                              0.1                  
  *N*                              20*D*                
  CMS                              *n*/4, if integer    
  *n*/2, if *n*/4 is not integer                        
                                                        
  CS                               *n*                  20*D* 
  *p*                              0.25                 
                                                        
  BBPSO                            *n*                  20*D*

###### 

Minimum NFE for the average best value to reach 1E-3, 1E-4, 1E-5, 1E-6, and 1E-7 from the known global minimum.

  Metaheuristic   *ε*          Phase equilibrium thermodynamic problem                                                                                      
  --------------- ------------ ----------------------------------------- ------------ ---------- ---------- ----------- ----------- ----------- ----------- ----------
  MAKHA           1E-3         *∞*                                       2164         601        121        301         481         4328        *∞*         10845
  1E*-*4          *∞*          45444                                     1803         121        602        4329        39493       *∞*         40006       
  1E*-*5          *∞*          *∞*                                       *∞*          **484**    903        *∞*         *∞*         *∞*         40488       
  1E*-*6          *∞*          *∞*                                       *∞*          2299       1806       *∞*         *∞*         *∞*         40729       
  1E*-*7          *∞*          *∞*                                       *∞*          **2783**   *∞*        *∞*         *∞*         *∞*         41211       
                                                                                                                                                            
  IFA             1E*-*3       *∞*                                       8820         1400       **40**     400         1600        19980       17300       **9440**
  1E*-*4          *∞*          24840                                     5200         **40**     3500       18880       36720       *∞*         **14480**   
  1E*-*5          *∞*          39600                                     *∞*          1040       9000       *∞*         52560       *∞*         **20080**   
  1E*-*6          *∞*          *∞*                                       *∞*          *∞*        15400      *∞*         73260       *∞*         25040       
  1E*-*7          *∞*          *∞*                                       *∞*          *∞*        *∞*        *∞*         *∞*         *∞*         30160       
                                                                                                                                                            
  CMAES           1E*-*3       *∞*                                       3961         **201**    41         **101**     **161**     3961        *∞*         *∞*
  1E*-*4          *∞*          5761                                      **201**      41         **101**    **2721**    5401        *∞*         *∞*         
  1E*-*5          *∞*          **7381**                                  **5801**     *∞*        **101**    **4801**    **7021**    *∞*         *∞*         
  1E*-*6          *∞*          *∞*                                       *∞*          *∞*        **101**    **8641**    **8821**    *∞*         *∞*         
  1E*-*7          *∞*          *∞*                                       *∞*          *∞*        **2401**   **11681**   **10621**   *∞*         *∞*         
                                                                                                                                                            
  ABC             1E*-*3       **9840**                                  8010         4100       60         850         2480        6570        **9750**    9560
  1E*-*4          *∞*          12150                                     9700         60         1950       7920        11790       **24959**   48236       
  1E*-*5          *∞*          17010                                     21500        1220       3850       34644       19530       *∞*         *∞*         
  1E*-*6          *∞*          *∞*                                       **58104**    12460      8650       *∞*         27990       *∞*         *∞*         
  1E*-*7          *∞*          *∞*                                       *∞*          33092      *∞*        *∞*         41850       *∞*         *∞*         
                                                                                                                                                            
  BA              1E*-*3       *∞*                                       **1980**     400        80         200         320         **1800**    *∞*         *∞*
  1E*-*4          *∞*          **2340**                                  400          80         200        5760        **3240**    *∞*         *∞*         
  1E*-*5          *∞*          *∞*                                       43200        6960       400        *∞*         16920       *∞*         *∞*         
  1E*-*6          *∞*          *∞*                                       *∞*          8400       3700       *∞*         *∞*         *∞*         *∞*         
  1E*-*7          *∞*          *∞*                                       *∞*          8440       *∞*        *∞*         *∞*         *∞*         *∞*         
                                                                                                                                                            
  MCSS            1E*-*3       87520                                     35640        800        **40**     400         640         26460       *∞*         21200
  1E*-*4          *∞*          80820                                     6200         **40**     1100       36800       104400      *∞*         34160       
  1E*-*5          *∞*          178920                                    78000        3520       3300       157600      *∞*         *∞*         64320       
  1E*-*6          *∞*          *∞*                                       *∞*          7880       11000      *∞*         *∞*         *∞*         *∞*         
  1E*-*7          *∞*          *∞*                                       *∞*          9120       49600      *∞*         *∞*         *∞*         *∞*         
                                                                                                                                                            
  CS              1E*-*3       41440                                     17460        7400       120        1700        7840        22140       16100       10640
  1E*-*4          **79200**    35100                                     11800        120        4100       26400       44820       30900       17680       
  1E*-*5          **109280**   61380                                     57400        1160       5500       52960       74700       **50300**   30000       
  1E*-*6          **132640**   *∞*                                       93400        3480       6700       93280       112500      **70300**   39760       
  1E*-*7          **155040**   *∞*                                       **131400**   4440       24100      140000      152820      **87300**   53840       
                                                                                                                                                            
  BBPSO           1E*-*3       *∞*                                       9900         2800       **40**     600         1760        7380        *∞*         **9440**
  1E*-*4          *∞*          14040                                     7400         **40**     1700       12640       11340       *∞*         17920       
  1E*-*5          *∞*          18720                                     21600        920        3300       *∞*         15300       *∞*         21440       
  1E*-*6          *∞*          *∞*                                       *∞*          **2200**   6000       *∞*         19080       *∞*         **23360**   
  1E*-*7          *∞*          *∞*                                       *∞*          5120       *∞*        *∞*         22500       *∞*         **28400**   

[^1]: Academic Editor: Xin-She Yang
