An Overview of the Application of Harmony Search for Chemical Engineering Optimization

Abstract

Harmony search algorithm and its variants have been used in several applications in medicine, telecommunications, computer science, and engineering. This article reviews the global and multi-objective optimization for chemical engineering using harmony search. The main features of the HS method and several of its popular variants and hybrid versions including their relevant algorithm characteristics are described and discussed. A variety of global and multi-objective optimization problems from chemical engineering and their resolution using HS-based methods are also included. These problems involve thermodynamic calculations (phase stability analysis, phase equilibrium calculations, parameter estimation, and azeotrope calculation), heat exchanger design, distillation simulation, life cycle analysis, and water distribution systems, among others. Remarks on future developments of HS and its related algorithms for global and multi-objective optimization in chemical engineering are also provided in this review. HS is a reliable and promising stochastic optimizer to resolve challenging global and multi-objective optimization problems for process systems engineering.

1. Introduction

Optimization problems are inherent for the process systems engineering of almost all units and operations associated to chemical engineering industry [1, 2]. Chemical engineers frequently face the scenario of identifying the best conditions to operate, for example, chemical reactors, mixing and separation equipment, polymerization systems, waste treatment and depollution processes, conventional and intensified separation sequences, mass-energy integration schemes, and supply chains. The resolution of these optimization problems allows to increase the economic profit, to reduce the energy consumption and its related environmental impacts, and to improve the sustainability and other target metrics/characteristics of the process under analysis. Therefore, the development and application of reliable numerical procedures for achieving the optimal process design and operation play an important role in the modern chemical engineering [2, 3].

Formally, the optimization refers to the implementation of mathematical and numerical procedures to calculate the best values of the decision variables (also known as optimization variables) that maximize or minimize one or more objectives of the analyzed system under the restrictions generated from a given scenario [4]. The optimization problems can be classified considering the characteristics of their objective function(s), decision variable(s), and constraint(s). For example, a problem is known as discrete, continuous, or combinatorial depending on the mathematical nature of variables that integrate the solution vector. The objective function properties determine if a given optimization problem can be classified as convex or non-convex, which can be also constrained by the restriction(s) of the search space via equality and inequality equations [5].

Chemical engineering optimization problems usually involve systems that are represented by non-linear complex (analytical and/or differential) equations derived from the mathematical models employed to describe the thermodynamics, transfer phenomena, and performance metrics. Several decision variables (discrete and continuous) that correspond to the degrees of freedom (i.e., independent variables) of the system can be involved in the problem solving. These optimization problems usually imply non-convex objective functions with several optimum values including the possibility of trivial and non-physical solutions. This context generates a continuous demand of reliable and effective optimization approaches to resolve challenging problems from chemical engineering applications.

Optimization techniques can be broadly classified as stochastic and deterministic [6–9]. Deterministic procedures like interval-based methods [10, 11], branch and bound [12, 13], primal-dual [14], and outer-approximation approaches [15] can guarantee theoretically the global optimum solution for some types of chemical engineering problems. However, they may require an extensive computational time to find the global solution (which is usually higher than those of stochastic optimization methods); besides, specific properties of objective functions and mathematical reformulations are needed to ensure the global optimum, thus limiting the resolution of complex problems [6, 16]. Alternatively, stochastic optimization methods can outperform several limitations of deterministic ones via the utilization of random search sequences coupled with heuristics and metaheuristics to generate a solution to the optimization problem [17, 18]. These optimizers are useful to resolve a variety of complex problems such as black-box algorithms, thus implying no restrictions in terms of the characteristics of objective function(s), problem dimensionality, and presence of constraints [19, 20]. A balance between diversification (i.e., exploration that is the ability of discovering promising areas to find the global optimum in the search space) and intensification (i.e., exploitation that utilizes the available information to guide the search process for improving the solution quality) is of great importance in these optimization algorithms [17, 21, 22]. A wide spectrum of stochastic optimization methods has been proposed and applied to solve chemical engineering problems [17, 23–26]. Some well-known techniques are cuckoo search (CS) [27], differential evolution (DE) [28], particle swarm optimization (PSO) [29], ant colony optimization (ACO) [30], tabu search (TS) [31], simulated annealing (SA) [32], genetic algorithm (GA) [33], and harmony search (HS) [34]. Novel strategies have been also introduced like grey wolf optimization (GWO) [35], flower pollination algorithm (FPA) [36], water cycle algorithm (WCA) [37], stochastic fractal search (SFS) [38], water wave optimization (WWO) [39], elephant herding optimization (EHO) [40], crow search algorithm (CSA) [41], sine cosine algorithm (SCA) [42], salp swarm algorithm (SSA) [43], butterfly optimization algorithm (BOA) [44], emperor penguins colony (EPC) [45], group teaching optimization algorithm (GTOA) [46], and smart flower optimization algorithm (SFOA) [47].

Since its publication, HS has been considered as a core reference method for researchers from different research and science fields. Several variants of this optimizer as well as its hybrid versions have been proposed and tested to solve a variety of global and multi-objective engineering problems [48, 49]. Results of these numerical studies have proved that key advantages of HS rely on its capability to handle discrete and continuous variables, its ability to escape from local minima, its independency of initial values for the decision variables, and its straightforward algorithm that utilizes simple operators and calculations [22, 48, 50]. Different reviews on HS have been published covering the following topics: HS algorithm and its most relevant variants, improvements, and hybridizations with other metaheuristics including applications in robotics, medicine, electrical engineering, image processing, data mining, civil engineering, power systems, manufacturing and design, scheduling, networking, agriculture, and water resource management [16, 48, 49, 51–54]. For example, Yoo et al. [55] published an overview about the application of HS in civil engineering. Askarzadeh [50] and Nazari-Heris et al. [56] carried out a literature review on HS-based methods to optimize electrical system problems. Yi et al. [57] presented a survey on HS and its applications on intelligent manufacturing. Zhang and Geem [58] provided a detailed study on the historical development of HS where the discussion focused on selected variants and hybridizations of this algorithm and their parameters, alternative initialization procedures, and constraint handling and multi-objective optimization techniques. Abualigah et al. [59] reviewed HS and its variants on clustering applications, while Nasir et al. [60] focused on the potential use of this stochastic optimization method for mathematics, computer science, and several engineering areas in three particular countries (i.e., China, South Korea, and Japan). Nasir et al. [61] presented a detailed overview about the combination and application of fuzzy logic theory within HS. However, the application of HS for global and multi-objective optimization in chemical engineering has not been covered in the available reviews despite its successful results in this engineering field.

The aim of this review is to provide the readers a broad perspective on the use of HS in a diverse set of optimization problems from the chemical engineering. Section 2 of this document briefly covers the ideas behind the development of HS, thus providing a comprehensive description of its algorithm and relevant variants. A discussion on HS algorithms for multi-objective optimization including a rich collection of algorithms is also provided in this section. Section 3 provides information on the application of HS to resolve optimization problems from chemical engineering such as the phase stability analysis and Gibbs free energy minimization of both reactive and non-reactive systems, the parameter identification of classical thermodynamic models like equations of state and local composition models, the calculation of azeotropes, critical points, and saturation conditions, the design of shell and tube heat exchangers and heat exchanger networks, the design and optimization of distillation sequences including intensified equipment (e.g., dividing wall column and thermally coupled reactive distillation), the life cycle analysis of energy generation systems, the design and optimization of water distribution systems, the multi-objective optimization of pump scheduling, the resolution of chemical equation balance, the parameter estimation in reaction kinetics, the chemical reactor design, the optimization and design of thermal cracking, and the resolution of optimal control problems. An overview of opportunities areas and current challenges to improve the performance of HS-based optimizers is given in Section 4. Finally, this review indicates that HS is a reliable and competitive algorithm to resolve a wide spectrum of global and multi-objective optimization problems from chemical engineering.

2. Fundaments of Harmony Search

2.1. Description of Basic Algorithm

In the musical improvisation, each musician plays his/her instrument spontaneously in a possible range of musical notes to achieve harmony together. When the musicians improvise, they could play a previously learned note, play something similar to that note and adjust it to the desired pitch, or play a new note. The quality of the improvised harmony is evaluated via an aesthetic standard, thus looking for quality and perfection in the performed musical piece. This music conception inspired Geem et al. [34] to develop the metaheuristic of HS. By comparing the musical improvisation with the optimization process, each decision variable (i.e., degree of freedom from the process under analysis) corresponds to a musician, the range of musical notes of an instrument refers to the range (i.e., search space) of the optimization variables, and, finally, the improvised harmony represents the solution vector for a given iteration [34] (see Figure 1). Just as the musical harmony is enhanced continuously, the improvement of the objective function value is expected at each performed iteration.

Figure 1 

Illustration of the analogy between musical improvisation and optimization concepts [62].

Mathematically speaking, an unconstrained single-objective (global) optimization problem can be defined assubject towhere is the objective function, is the solution vector, is a decision variable subject to lower and upper bounds and , respectively, and is the problem dimension (i.e., number of decision variables). For the case of constrained problems, the optimization of objective function given by equation (1) must satisfy the following restrictions:where q and p are the number of equality and inequality constraints, respectively.

Herein, it is convenient to note that constraints given by equations (3) and (4) can be easily handled by HS using a penalty approach [63, 64]. Constrained and unconstrained global optimization problems are commonly formulated during the chemical engineering practice and they include the parameter identification of thermodynamics, kinetic and transfer phenomena models [1, 65–67], phase equilibrium calculations [68, 69], and design of purification systems using a single objective or target metric [70–72], among others. HS is an effective optimizer to resolve a wide spectrum of global optimization problems from chemical engineering [49, 73]. This method has four parameters to control the search process to minimize or maximize the tested objective function [34]: pitch bandwidth (bw), pitch adjusting rate (PAR), harmony memory considering rate (HMCR), and harmony memory size (HMS). HS optimizer comprises four main steps: (1) initialize the harmony memory (HM), (2) improvise a new harmony, (3) update HM, and (4) perform the iterative cycle until the stopping condition is satisfied. A brief description of these steps is provided below.

Step 1. Initialize HM. In this stage, each decision variable of the solution vectors is randomly generated and stored in HM with their respective objective function valuesEach candidate solution vector is evaluated and sorted in HM according to its value of objective function: .

Step 2. Improvise a new harmony. In this stage, a new harmony (i.e., new solution vector) is obtained using three numerical operators: random selection, pitch adjustment, and memory consideration. This improvisation process is mainly controlled by two probabilistic parameters, PAR and HMCR, where 0 ≤ HMCR, PAR ≤ 1, that are applied sequentially in the algorithm to generate a new set of candidate solutions. First, the new values of decision variables are obtained according towhere rand is a random number within the variable bounds. In the memory consideration mechanism, HMCR is the probability of selecting one element from HM as the new decision variable value. The random selection occurs with probability (1 – HMCR) and the new decision variable value is calculated randomly. If is obtained via the memory consideration, the pitch adjustment procedure is performed according to the following operator:where rand₀ is a random number between 0 and 1. The improvisation step considers all the harmonies stored in HM.

Step 3. Update HM. The objective function is evaluated using ; if its value is better than the worst stored in HM, , the new solution vector will substitute its place in HM.

Step 4. Iterative cycle. Steps 2 and 3 are repeated until the specified termination criterion is satisfied.
Illustratively, Figure 2 presents a simplified diagram of the HS algorithm. It is convenient to recall that the performance of a stochastic optimization algorithm relies on the trade-off between exploration and exploitation stages [74, 75]. HMCR and PAR affect the success rate during the global optimization. PAR parameter controls the local search where a low value with a narrow bw can delay the algorithm convergence since the exploration is limited to a small region of the search space. In contrast, a high PAR value with a large bw enables HS to escape from areas potentially close to the local optimal solutions [53]. The exploitation stage is mainly governed by HMCR, which has been considered as a key parameter that determines the level of elitism. A low value of this parameter implies a slow convergence to the optimal solution. On the contrary, a high HMCR value favors the algorithm convergence but at the same time generates the possibility that HS could be trapped in a local solution [53]. Both parameters influence the performance of HS to find global or local solutions, and an adequate balance between exploitation and exploration stages should be obtained via their tuning [21, 22]. The determination of the best HS algorithm parameters has proved to be problem dependent. Therefore, previous studies and preliminary calculations are useful as a starting point to define these values for solving the optimization problem at hand. Some authors have indicated that HS effectiveness could not be very sensitive to the values of these parameters for some optimization problems and, under these circumstances, it is not necessary to perform an exhaustive adjustment of them [76, 77].

Figure 2 

Flowchart of harmony search algorithm.

2.2. Variants of Harmony Search

Since the initial version of HS algorithm was introduced [34], a large number of variants have been developed by different researchers with the aim of improving its performance. From a general perspective, the modifications of this optimization method are usually made in the algorithm structure or in the way that its parameters are determined. Also, a wide variety of hybrid HS-based approaches can be found in the literature. A brief description of relevant HS variants is provided below.

2.2.1. Improved Harmony Search

Mahdavi et al. [78] developed the improved harmony search (IHS). The main difference between HS and IHS relies on the improvisation step. HS uses fixed PAR and bw parameters, while IHS calculates them dynamically according to the iteration number viawhere subscripts min and max indicate the minimum and maximum values of bw and PAR parameters and iter is the actual iteration. PAR increases linearly with the maximum number of iterations (iter_max), while bw decreases exponentially.

2.2.2. Global-Best Harmony Search

The global-best harmony search (GHS) was introduced by Omran and Mahdavi [79]. This method incorporates the concept of swarm intelligence used in PSO. A swarm is a set of particles in PSO where each solution vector is presented by a particle. The location of the vector in the search space is given by its best fitness and the best global solution found by all the swarm particles [29]. In GHS, the improvisation step was modified implying that will take the value of the corresponding element of the best harmony in HM in case of pitch adjustment. The bw parameter is omitted and PAR is calculated according to equation (8). A comparison of GHS and HS improvisation procedures is provided below (see Table 1).

Table 1 

Comparison of the improvisation step in HS and GHS.

2.2.3. Novel Global Harmony Search

A new HS algorithm called the novel global harmony search (NGHS) was reported by Zou et al. [80]. Based on swarm intelligence, NGHS adds two relevant features to the improvisation step: genetic mutation and position updating. The original HS parameters (i.e., HMCR and PAR) are substituted by a genetic mutation probability to avoid local solutions caused by premature convergence. Also, replaces even if it has worst fitness when updating HM. The improvisation step of NGHS is described below (see Table 2).

Table 2 

Improvisation step of the NGHS algorithm.

2.2.4. Self-Adaptive Global-Best Harmony Search

The self-adaptive global-best harmony search (SGHS) was based on GHS [81] and incorporates a dynamic parameter adjustment mechanism to define bw, PAR, and HMCR. Since defining these algorithm parameters is problem dependent, its adaptative adjustment makes SGHS suitable to solve optimization problems with different characteristics. Initially, PAR and HMCR are calculated based on a fixed mean value and standard deviation assuming a normal distribution. After certain number of iterations, these mean values are recalculated and HMCR and PAR are updated. Regarding bw, its value is dynamically reduced as the iterations increase:

2.2.5. Improved Global-Best Harmony Search

El-Abd [82] introduced an enhanced version of GHS, which was called the improved global-best harmony search (IGHS). In the improvement step, the memory consideration is performed as follows:where a Gaussian distribution with mean 0 and standard deviation of 1 is used to calculate the random number . On the other hand, the pitch adjustment mechanism is given bywhere  [−1, 1]. In this algorithm, PAR decreases linearly over iterations, while HMCR remains constant during the optimization process.

2.2.6. Improved Differential Harmony Search

Wang et al. [83] developed the improved differential harmony search (IDHS) algorithm. IDHS integrates two solution vector generation strategies commonly used in DE: DE/best/1/bin and DE/rand/1/bin. The former is incorporated in memory consideration and the second is applied to random selection. bw is removed in the pitch adjustment to include a swarm intelligence concept where takes the value of as in GHS. HMCR and PAR are dynamically adapted according to

2.2.7. Enhanced Self-Adaptive Global-Best Harmony Search

Luo et al. [84] developed the enhanced self-adaptative global-best harmony search method (ESGHS). A parameter-free scheme was proposed where HMCR can be updated dynamically according to a normally distributed random number in a problem dimension dependent range. PAR is calculated using the current and maximum number of iterations and bw is defined bywhere superscripts h and m refer to random integers within [1, HMS]. The pitch adjustment mechanism of ESGHS is given by

If at any point, the following equation is used for the pitch adjustment:

In the random selection (1 − HMCR), the new harmony element is obtained as follows:where a Gaussian distribution with mean and standard deviation is used to calculate the random number rand_G.

2.2.8. Selective Acceptance Novel Global Harmony Search

The selective acceptance novel global harmony search (SANGHS) was proposed by Li et al. in [85] which is among the latest published articles regarding the development of HS variants. In this algorithm, an additional selective acceptance mechanism is incorporated besides the two mechanisms involved in NGHS (i.e., position updating and genetic mutation), which is inspired by the SA concept. An acceptance probability is calculated at each iteration to decide if the new solution is accepted or rejected following the next criterion:

2.2.9. Other Recent HS Algorithm Developments

The improved differential-based harmony search algorithm with linear dynamic domain (ID-HS-LDD) was introduced by Zhu et al. [86]. Two novel strategies are used in this algorithm: an improved differential-based technique (DE/best/2/bin) that replaces the original bw parameter and a dynamic change in the search space to avoid stagnation in local optima regions. These authors also introduced an additional parameter that modifies the scheme to generate new harmonies with the aim of enhancing the convergence speed.

Zhao et al. [87] developed a new adaptive HS algorithm (aHSDE) based on periodic learning, differential mutation, and linear population size reduction. This algorithm considers all the solutions stored in HM, including the best solution, to generate the new harmony elements using a differential mutation operator (DE/best/2). In order to balance the convergence rate and diversity, a linear reduction strategy is applied to handle HMS in terms of the maximum and minimum values of HMS and iterations. Also, the concept of period learning was used to dynamically adjust PAR and F (a scaling factor in the differential mutation operator) with the objective of improving the algorithm adaptability.

Jeong et al. [88] presented an advanced-parameter-setting-free strategy to handle HS parameter definition. HMCR parameter was calculated with respect to the maximum and current numbers of iterations, as well as the number of optimization variables. A sigmoid function was involved in HMCR calculation to maintain its value within 0.5 and 1. PAR parameter was calculated with respect to HMCR, and the sigmoid function was also applied to a dimensionality dependent term.

The most relevant features of HS variants reviewed in this paper are summarized in Tables 3–5. In particular, Table 3 contains those algorithms based on the HS parameter modifications, and Table 4 includes algorithms based on the HS structure modification, while Table 5 provides information on HS hybridizations. Figures 3 and 4 provide the grouping of different HS-derived algorithms according to the changes in algorithm structure, parameter determination, hybrid techniques, and those techniques to address optimization problems with discrete variables.

Table 3 

Harmony search variants based on the operators to calculate the algorithm parameters.

Table 4 

Harmony search variants based on the algorithm structure modification.

Table 5 

Harmony search variants based on hybridization with other metaheuristics.

Figure 3 

Harmony search variants grouped by the structure-based modifications.

Figure 4 

Harmony search variants grouped by the parameter-based modifications, hybrid approaches, and multi-objective methods.

2.2.10. Multi-Objective HS Algorithms

Several chemical engineering problems can involve several objectives that are usually in conflict and should be optimized simultaneously [1, 190]. The optimal design of conventional and intensified separation sequences to maximize the purity of target compound(s) and to minimize the energy consumption [191, 192], the integration of mass and energy in different processes [193–195], and the simultaneous minimization of economic (i.e., capital and operating costs) and environmental objectives [196] are examples of multi-objective optimization problems in the context of chemical engineering. The resolution of these problems can be performed using multi-objective HS algorithms. Formally, an optimization problem with two or more objectives can be defined aswhere m is number of objective functions to optimize simultaneously. This multi-objective problem can be also subject to inequality and equality constraints given by equations (2)–(4). The determination of Pareto front(s) is the main characteristic of the resolution of multi-objective optimization problems, which are used to analyze the relationships between conflicting objective functions and the corresponding values of decision variables [197].

Stochastic multi-objective optimization methods have been of great interest in process systems engineering due to their advantages to handle both linear and non-linear equality and inequality constraints and their capabilities to escape from local solutions and to employ any type of objective functions and decision variables [198].

The first applications of HS to solve multi-objective optimization problems were carried out by Geem [199] and Gemm and Hwangbo [200] for the scholar bus routing problem and satellite heat pipe design, respectively. Afterwards, Sivasubramani and Swarup [201] developed the multi-objective harmony search (MOHS) algorithm. In this method, the PAR and bw parameters change at each generation according to a fixed maximum and minimum values and the current and maximum number of iterations. To update HM, a non-dominated sorting and ranking mechanism and a dynamic crowding distance strategy are used [202, 203]. The fuzzy membership approach is employed to find a balanced Pareto front solution between the conflicting objectives.

Ricart et al. [204] presented two algorithms to solve optimization problems with two or more objectives: MOHS1 and MOHS2. The first method employs a ranking assignment method proposed by Fonseca and Fleming [205] where the best trade-off solutions are stored in HM by considering their rank. MOHS2 algorithm utilizes an alternative memory at each iteration, which is combined with the original memory, to select a limited number of solutions for the next iteration. Both algorithms were compared with the non-dominated sorting genetic algorithm II (NSGA-II) on ZDT functions.

Pavelski et al. [206] proposed a series of multi-objective optimization algorithms based on HS, IHS [78], GHS, and SGHS. These methods were based on the NSGA-II concept, so the novel algorithms were called non-dominated sorting harmony search (NSHS), NSIHS, NSGHS, and NSSGHS, respectively. For all the methods, the quantity of solution vectors in the memory is duplicated, and the non-dominated sorting with elitism and crowding distance is applied to update HM. The performance of these algorithms was assessed according to the Friedman test via the solving of the Congress on Evolutionary Computation 2009 multi-objective benchmark set.

Nekooei et al. [207] developed an improved version of NGHS called improved multi-objective harmony search (IMOHS). This multi-objective method was based on the combination of all the objective functions involved in the problem into a single objective through a weighted sum following the approach of previous works [208, 209]. As in MOHS, the search process is governed by non-dominated ranking and sorting. An additional memory known as archive is implemented to keep a record of non-dominated harmonies at each iteration; meanwhile, the dominated harmonies are stored in HM. An update process based on a specified rank of dominated harmonies is considered to maintain an adequate relationship between exploitation and exploration stages of the algorithm.

Sabarinath et al. [210] proposed the multi-objective parameter adaptative harmony search (PAHS). PAHS algorithm was adapted to handle multi-objective optimization problems using the weighted sum approach. On the other hand, Prajapati and Chhabra [211] developed the many-objective discrete HS (MaDHS) algorithm. They proposed to solve discrete optimization problems with more than three objectives based on non-dominated solution ranking instead of the Pareto approach. Molina-Pérez et al. [212] developed multi-objective HS (MOHSg). This method includes a modification in the pitch adjustment step using a crowding distance by genotype to improve the exploration/exploitation balance. Hesar et al. [213] developed a quantum multi-objective harmony search algorithm (QMOHSA) based on HS, evolutionary concepts (i.e., non-dominated sorting and elitism), and quantum computing. The exploration/exploitation balance of this algorithm is handled using quantum computing concepts such as qubits and quantum gates and HS operators (i.e., PAR and bw). This multi-objective algorithm was compared with other three multi-objective techniques on several benchmark functions and results showed its better performance. Figure 4 summarizes the multi-objective methods based on HS metaheuristic.

3. Applications of Harmony Search in Chemical Engineering Optimization

The process systems engineering deals with chemical engineering problems that can be resolved using global and multi-objective optimization approaches. HS has been widely utilized in this field, and this section aims to provide a broad perspective on its chemical engineering applications. For illustration, Figure 5 shows the number of articles found per year on the development of HS variants as well as on its application in chemical engineering. A concise description of the optimization problem, involved variables, problem constraints, and relevant remarks of representative applications is provided in the following sections where the summarized information is given in Tables 6–12 for an easy reference for interested readers.

Figure 5 

Number of articles published per year on harmony search variants and their applications in chemical engineering optimization.

Table 6 

Harmony search applications in chemical engineering optimization.

Table 7 

Problem description of thermodynamic calculations.

Table 8 

Problem description of the design of heat exchangers.

Table 9 

Problem description of distillation design and optimization.

Table 10 

Problem description of life cycle analysis.

Table 11 

Problem description of water distribution system design.

Table 12 

Problem description of other applications in chemical engineering.

3.1. Thermodynamic Calculations for Process Systems Engineering

Several thermodynamic calculations for the process systems engineering of reactors, separation and purification units, and processes can be formulated as optimization problems. They include the parameter estimation of thermodynamic models, phase equilibrium calculations and stability analysis in non-reactive and reactive systems, and prediction of azeotropes, saturation conditions, and critical points [69, 253]. A brief description of those thermodynamic problems solved with HS-based methods is provided below.

3.1.1. Phase Stability Analysis

This analysis allows to identify if the thermodynamic state of a mixture is stable at given operating conditions. A given mixture at certain feed composition z, pressure, and temperature is stable if the Gibbs free energy surface lies above the tangent plane to such surface at z for all other composition points [254–256]. This calculation is a fundamental stage to be performed prior to the phase equilibrium calculations of non-reactive and reactive multi-component mixtures. The tangent plane distance function (TPDF) is the objective function to be globally minimized to resolve this thermodynamic problem [255]. This optimization problem is constrained by the material balance restrictions and non-negative bounds on the compositions that are the corresponding decision variables (see Table 7). The identification of the global minimum of TPDF is the correct solution of this thermodynamic problem to avoid improper conclusions on the phase equilibria behavior of analyzed mixtures, which directly affect the process design of separations systems and multi-phase units.

3.1.2. Phase Equilibrium Calculations

They are relevant for the modeling of separation processes (e.g., distillation and extraction) and multi-phase reactors. The objective in phase equilibrium calculations is to establish the correct type, number, and composition of phases present at thermodynamic equilibrium. The Gibbs free energy minimization is the optimization problem that is resolved for non-reactive and reactive systems [69]. It is a constrained global optimization problem due to the mass balance restrictions where the decision variables are the number of phases and their compositions. The chemical reactions increase the dimensionality and complexity of the optimization problem for reactive mixtures where additional constraints related to the chemical equilibrium should be considered [257]. The location of the global minimum of Gibbs free energy function is essential for phase equilibrium calculations because it corresponds to the only correct and desirable solution for performing reliable process design and modeling [256]. Dohrn and Pfohl [258] illustrated the significant impact that small uncertainties in phase equilibrium calculations could have on the distillation column design. For instance, 5% underestimation of the separation factor generated 100% overestimation of the column height, thus significantly increasing the equipment cost and energy requirements. Hence, the phase equilibrium calculations must be performed efficiently and reliably to avoid process design errors.

3.1.3. Parameter Identification of Thermodynamic Models

In general, the thermodynamic models have adjustable parameters that must be obtained in order to perform the calculation and estimation of properties of pure components and mixtures [67, 259]. These model parameters can be estimated from the minimization of an objective function that usually is defined in terms of available experimental information (e.g., equilibrium compositions and saturation pressures) for the system (e.g., pure component or mixture) under analysis. Different objective functions can be formulated for this optimization task where the error-in-variable and the classical least-squares formulations are the common ones [224]. These objective functions are usually non-convex where unconstrained and constrained global optimization problems should be resolved.

3.1.4. Estimation of Azeotropes, Critical Points, and Saturation Conditions

Special phase equilibrium calculations are also required for the process systems engineering of separation systems. They include the calculation of homogeneous and heterogeneous azeotropes [260, 261], critical transitions [262], and bubble and dew points [263, 264]. For instance, the azeotropes occur in a boiling mixture with one or two liquid phases when the composition of the vapor is the same as the overall composition of the liquid phase(s). When one liquid phase is at equilibrium with a vapor phase, the azeotrope is called homogeneous azeotrope. In case of two liquid phases at equilibrium with a vapor phase, the azeotrope is known as heterogeneous [261]. These special phase equilibrium cases can also occur in systems where chemical reactions take place and they are called as reactive azeotropes [260]. The azeotropy is relevant for the operation of separation systems based on vapor-liquid equilibrium and, consequently, it can occur in many industrial applications. One well-known example is the ethanol + water distillation process. The homogeneous azeotrope of this mixture imposes a limit to the purity degree that can be reached via distillation. The capability to predict this condition reliably in multi-component mixtures is essential for the separation process design [260]. On the other hand, the calculation of critical transitions of multi-component mixtures is also relevant for the operation and design of separation systems [262], while the bubble and dew point calculations are required for the analysis of phase diagrams of mixtures [265]. All these calculations are usually performed by solving a set of non-linear equations that define the thermodynamic equilibrium condition to be satisfied, which usually corresponds to the equality of chemical potentials [263, 266]. Therefore, a global optimization approach can be used to resolve this set of non-linear equations. For example, the optimization variables for azeotrope calculations are the azeotrope composition and its temperature or pressure. However, it should be noted that these thermodynamic calculations are characterized by the presence of several solutions (stable and unstable from the perspective of phase stability analysis) or even no solutions for the problem at hand [267]. Interested readers can consult the studies of Henderson et al. [268]; Bonilla-Petriciolet et al. [263]; Bonilla-Petriciolet et al. [269]; and Nakazawa et al. [270] for more details on the formulation of these optimization problems.

HS-based algorithms have been applied to resolve some of these thermodynamic calculations. In particular, Bonilla-Petriciolet et al. [218] carried out the parameter estimation of local composition models for calculating VLE in binary systems. HS was used to minimize an objective function defined in terms of activity coefficients following the least-squares formulation. The activity coefficients were calculated using several local composition models: Wilson, UNIQUAC, and NRTL. The global optimization of these parameter identification problems was performed using the following HS algorithm parameters: PAR = 0.75, HMCR = 0.25, and HMS = , where n_var was the number of decision variables. The decision variables were the binary interaction parameters of Wilson, UNIQUAC, and NRTL models. A quasi-Newton method was incorporated at the end of HS optimization to improve the final solution accuracy. According to several performance metrics, these authors concluded that HS was a promising metaheuristic for the VLE modeling.

In other study, Bonilla-Petriciolet et al. [220] compared the performance of GHS, IHS [78], and HS in the phase equilibrium modeling of systems without chemical reactions. Phase equilibrium calculations and stability analysis were performed via the global minimization of the Gibbs free energy and TPDF, respectively. Both calculations were resolved as unconstrained optimization to reduce the problem dimension [271] (see Table 7). Nine mixtures from two to ten components showing LLE and VLE were analyzed using NRTL and SRK models. The following parameters were used for HS-based algorithms: bw = UB_i − LB_i, PAR = 0.75, and HMCR = 0.5 for HS and GHS and bw_min = 0.001, bw_max = UB_i − LB_i, PAR_min = 0.5, PAR_max = 0.95, and HMCR = 0.5 for IHS. For all cases, HMS was defined as . Two termination criteria were used: a maximum number of iterations without improvement and iter_max. A quasi-Newton method was incorporated to GHS and IHS to assess their performance. A detailed analysis based on performance profiles [272] was provided for these HS algorithms. In general, IHS and GHS were more reliable than HS to resolve these thermodynamic calculations.

Merzougui et al. [222] compared the performance of HS and GA to estimate NRTL binary interaction parameters. Twenty-one ternary LLE mixtures were studied. The parameter estimation problem was formulated via the minimization of the squared error sum between calculated and experimental compositions of all components in both liquid phases. A local search based on a quasi-Newton method was implemented to improve the accuracy of the final solution found by HS. The parameters of HS were determined by a sensitivity analysis and they were defined as HMCR = 0.9, PAR = 0.4, HMS = 10, and iter_max = 16,000. The authors concluded that the binary interaction parameters found using HS will provide more accurate phase equilibrium results than those obtained using GA.

Bonilla-Petriciolet et al. [223] assessed the performance of GHS, IHS [78], and HS algorithms in the solution of phase equilibrium problems of reactive systems using the Gibbs free energy minimization approach. The constraints of this optimization problem were handled using a formulation that reduces the number of decision variables and numerical effort. The penalty function approach was utilized to handle these constraints. Eight reactive mixtures with LLE and VLE were studied using several thermodynamic models (i.e., NRTL, Wilson, ideal gas, UNIQUAC, and Margules). The parameters of all optimization methods were set according to results from previous studies [220]. IHS and GHS were able to escape from the infeasible zones easier than HS, and they also showed better exploitation characteristics. These authors concluded that the performance of HS and its variants was problem dependent, but they were reliable to solve this type of thermodynamic problems.

Bonilla-Petriciolet [224] employed HS to perform the parameter estimation of UNIQUAC, Wilson, and NRTL models for VLE modeling of binary mixtures. Isothermal and isobaric data were used to formulate the error-in-variable and least-squares objective functions. The least-squares formulation was based on activity coefficients, and the decision variables corresponded to the binary interaction parameters of NRTL, UNIQUAC, and Wilson models. On the other hand, the Wilson model was utilized in the error-in-variable formulation where the decision variables corresponded to the true values of the liquid composition and temperature to carry out the corresponding data reconciliation. bw_i = UB_i − LB_i, PAR = 0.75, HMCR = 0.25, and HMS =  were the parameters used for HS. This author concluded that HS showed a better performance in comparison with other stochastic algorithms like GA and PSO in terms of its success rate for the phase equilibrium data correlation using the least-squares formulation. However, the reliability of HS was reduced in the parameter estimation with the error-in-variable objective function due to its problem complexity. Therefore, a hybrid HS variant (HHS) incorporating SA and DE operators was introduced and compared with HS, GHS, and IHS [78], thus improving the performance of original HS.

Other studies have published the application of HS on the modeling of ionic liquid properties, which are considered as outstanding solvents for extraction operations and other processes [273, 274]. Ionic liquids (ILs) are electrolyte solutions with interesting physicochemical properties like very low vapor pressures, high conductivity, and thermal stability [275]. Bonilla-Petriciolet et al. [225] carried out the evaluation and comparison of several stochastic global optimization techniques including HS to model the mean activity coefficients of different ammonium electrolytes in aqueous solution using the eNRTL model proposed by Chen et al. [276]. The modeling of activity coefficients of electrolyte solutions using this local composition model is a challenging global optimization problem [277]. A least-squares formulation in terms of activity coefficients was defined to determine the binary interaction parameters of these thermodynamic model via global optimization. The performance of HS was assessed using performance profiles, success performance, success rate, and global success ratio. Results with and without a local optimizer were presented. HS algorithm was implemented with HMCR = 0.25, PAR = 0.75, and different stopping conditions. The parameter estimation process for this application was found to be challenging due to the presence of comparable optimums, local minimums, and saddle points. HS showed a success rate from 0 to 100% for finding the global optimum in this application.

Thermodynamic systems containing biofuels can show different types of phase equilibrium (e.g., LLE, VLE, and VLLE) due to the chemical nature of their components. The parameter estimation problem for this type of systems is a challenging global optimization problem [231]. Therefore, a comparative study on the performance of HS, GA, and backtracking search optimization algorithm (BSOA) for the binary interaction parameter identification of NRTL and UNIQUAC models using thirty LLE ternary systems relevant in biodiesel production was carried out by Merzougui et al. [231]. The objective function was defined in terms of calculated and experimental compositions of these LLE systems. In other study, Merzougui et al. [233] analyzed the numerical performance of HS, SA, flower pollination algorithm (FPA), a FPA variant, GA, and BSOA in the parameter estimation of NRTL and UNIQUAC binary interaction parameters in several relevant mixtures from the food industry. The global optimization was performed via the minimization of an objective function defined in terms of LLE concentrations of ternary and quaternary mixtures.

Fernández-Vargas et al. [234] studied the effect of six termination criteria in the performance analysis of different stochastic methods, including HS, in the global minimization of thermodynamic problems. Thirty-two problems from local composition model parameter estimation, Gibbs free energy minimization with and without chemical reactions, and phase stability analysis were studied. An unconstrained approach for the global minimization of TPDF and the Gibbs free energy for non-reactive systems with LLE and VLE was used. On the other hand, a constrained approach for the Gibbs free energy minimization with chemical equilibrium was applied, while the parameter estimation problems were formulated using least-squares for VLE data. The parameters of HS for resolving these problems were as follows: HMS = , HMCR = 0.5, PAR = 0.5, and iter_max = . Results showed that the success rate of each algorithm can change significantly according to the termination criterion utilized in these thermodynamic calculations. The error-based stopping criterion showed the best results followed by the improvement-based termination. Gibbs free energy minimization with chemical equilibrium and the parameter estimation problem were the most challenging calculations for tested stochastic optimization methods. These authors recommended the study of other termination conditions to improve the performance of stochastic optimizers in thermodynamic calculations, particularly in the parameter estimation problem.

GWO, WCA, and a combined algorithm from HS and GWO were employed for the global minimization of TPDF [246]. Several mixtures with LLE and VLE modeled using SRK, UNIQUAC, and NRTL were employed in this study. The incorporation of PAR parameter improved the exploration stage of GWO, thus allowing the identification of promising areas to find the TPDF global optimum. This study proved that the numerical operators of HS can be incorporated in the algorithm of other metaheuristics to enhance their numerical performance in thermodynamic calculations.

Supercritical fluid extraction is used as an alternative attractive approach for the decaffeination of coffee and tea leaves and extraction of essential oils and flavors from species and herbs, among other industrial applications [278]. A fluid is in a supercritical state when its temperature and pressure are above its critical value and the generated fluid shows interesting thermodynamic properties [279]. An adequate understanding of phase equilibrium and solubility data of supercritical fluids is essential for extraction processes [280]. Amaya et al. [243] performed the parameter estimation of Peng–Robinson (PR), Van Laar, NRTL, and UNIQUAC models for the VLE of supercritical-CO₂ + α-pinene system using the SFHS algorithm. This global optimization problem was resolved using the following SFHS parameters: HMS = 5, PAR = 0.8, HMCR = 0.9, bw_ini = 0.5, bw_min = 1000, bw_max = 2, bw_sat = 1000, and amplitude constant C_bw = 1. Results showed that SFHS was a promising technique for phase equilibrium modeling with supercritical fluids.

Dew point calculations in binary systems with double retrograde vaporization are characterized by the presence of three or four dew point values for a single composition [281]. This phenomenon only occurs in conditions near the critical temperature of the system and its calculation is challenging. This problem is particularly difficult to solve because several local optimums exist including trivial solutions with no physical meaning [282]. The system of non-linear equations derived from the isofugacity criterion was integrated in an objective function to perform the global optimization where the decision variables were the composition and pressure. Platt et al. [245] employed HS and DE to calculate dew point pressures in binary double retrograde vaporization systems and reactive azeotropes. The authors mentioned that solving a non-linear system of equations via an optimization approach is also challenging. The control parameters of both metaheuristics were tuned to achieve reliable results. This parameter tuning is relevant because the numerical performance of each algorithm is problem dependent and it should be performed in detail for some challenging optimization problems [69].

Reggab et al. [247] regressed NRTL and UNIQUAC binary interaction parameters for the modeling of hexane + 1-propanol + water mixture using HS and other metaheuristics. The objective function was defined in terms of mass fractions of every component on each phase. The performance of these metaheuristics was compared with and without a local solver (i.e., quasi-Newton or Nelder–Mead). The incorporation of a local optimization method reduced the average root mean squared deviation by 20%.

Smejkal et al. [248] carried out the performance analysis of several metaheuristics (HS, DE, CS, covariance matrix adaptation, and elephant herding optimization) and the Newton–Raphson method in phase stability analysis. TPDF in terms of the Helmholtz free energy density was minimized as an unconstrained problem. TPDF minimization was done with the following HS parameters: iter_max = 5E10⁴, bw = 0.001, PAR = 0.5, HMCR = 0.85, and HMS = 6. Recently, Fonseca-Pérez et al. [249] proposed a benchmark set of phase stability problems to properly evaluate the performance of current and new metaheuristics. A classification of phase stability problems was proposed via a numerical analysis on the performance of several stochastic optimization methods including HS to globally minimize TPDF. This minimization problem was solved as an unconstraint problem employing HS with the following parameters: PAR = 0.75, HMCR = 0.5, and HMS = . Global performance results placed HS in the third position of tested stochastic optimizers.

3.2. Design of Heat Exchangers

Heat exchangers refer to equipment that allows heat transfer between two or more fluids at different temperatures where different equipment can be used [283–286]. For example, the plate-fin heat exchangers (PFHEs) are compact instruments mainly composed of a set of stacked flat plates with corrugated fins [287]. The heat exchange takes place through the transition of a fluid between the corrugated fins. PFHE design involves an integrated analysis based on thermodynamics, mechanical design, transport phenomena, and cost estimation. Another example of this design problem is related to the shell and tube heat exchangers (STHEs) (see Figure 6). This equipment is commonly utilized in industrial processes due to its variety of architectures for heat transfer operations involving steam generators, evaporators, condensers, or water boilers. STHEs are able to operate at high temperature conditions and pressures. Its large heat transfer surface increases the equipment efficiency; besides, its installation, cleaning, and maintenance are relatively easy to perform [289]. STHE design also involves several variables that are defined considering the physical properties and geometric configuration of the exchanger. The main objective of exchanger design is to obtain a heat transfer coefficient that contributes to an efficient operation besides optimizing other attributes of the equipment operation. This analysis entails a variety of possible objective functions to be optimized where the minimization of the capital cost is the most common. Several studies have defined different design targets (e.g., the weight of the heat exchanger) for the heat exchanger design that can be handled via a multi-objective approach such as maximizing the effectiveness and simultaneously minimizing the annual cost [290, 291]. Overall, the optimum design of this type of equipment implies to minimize its operating cost, which involves the exchanger energy consumption and the capital cost including the design, materials, installation, testing, manufacturing, and shipping costs. This design problem has usually equality and inequality constraints that increase the complexity of problem resolution and equipment design [292]. Hence, reliable numerical tools are required to obtain an optimal design for heat exchangers where HS has proved to be useful.

Figure 6 

Illustration of shell and tube heat exchanger diagram [288].

Fesanghary et al. [252] compared the STHE design with HS and GA via the minimization of an objective function based on the capital and operating costs of the equipment for a required thermal load. This case of study considered the cooling of an oil through a water flow. First, a global sensitivity analysis was performed to identify the parameters with the highest impact on the STHE cost. It was determined that the inside shell and outside tube diameters were the main parameters. The number of sealing strips was discarded because it had the minor impact on the exchanger cost. Therefore, the objective function was formulated considering nine discrete and continuous variables (e.g., inside shell and outside tube diameters, type of material, and baffle cut). The parameters used for HS were as follows: bw ∈ (0.01 – 4), PAR ∈ (0.2 – 0.85), HMCR = 0.6, and HMS = 10. The evaluation of several configurations of the exchanger geometry was carried out. HS and GA found solutions near to the optimum value where the difference between the global optimum and the solution found by HS was 0.2%, while GA obtained a 0.7% deviation.

Khorasany and Fesanghary [217] developed a two-level hybrid approach using HS and SQP to minimize the total annual cost of a heat exchanger network. The objective function involved several terms associated to the utility and capital costs. The optimization variables were defined as the heat load of each exchange unit and the stream split fractions. Several constraints were considered in this problem: the overall heat balance for every stream, minimum temperature change, mole fraction constraints, and non-negative bounds on the optimization variables. The penalty function method with a static penalization was used to manage these constraints. The first stage of this hybrid approach was to determine the structure of the heat exchanger network (i.e., the number of heat exchangers and branches in each stream). Discrete variables were involved in this stage, and HS was used to handle them. In the second stage, the stream split fractions and the heat load of the exchangers (which were continuous variables) were optimized using HS with SQP to obtain the minimum total annual cost. A probability parameter was defined to decide whether to perform or not the SQP-based local search during the improvisation step using the new harmony as starting point. Once the termination criterion was satisfied, local optimization was performed to improve the solution quality. This algorithm was tested on benchmark and real-world problems. These calculations showed stable convergence and high efficiency of this algorithm.

Doodman et al. [215] evaluated the IHS [78] performance to minimize the total design cost of air-cooled heat exchangers. A sensitivity analysis was carried out to identify the influence of design parameters on the equipment cost and to reduce the optimization problem dimension. The fin number, fan type, and tube pattern were eliminated from the optimization process since they did not show a significant influence on the exchanger cost. The optimal geometry of the exchanger was described by the air velocity, tube outer diameter, tube material and length, number of tube passes, number of tubes per row, and fin height. The performance of IHS was compared with GA, thus finding similar results. The main differences between these optimizers were observed in the number of tubes passes and tubes per row, as well as in the tube length (i.e., 2 and 3 tube passes with 24 and 20 tubes per row for GA and HS, respectively). However, the solution precision of HS was better than GA.

Yousefi et al. [226] proposed an improved HS algorithm to optimize PFHE design. IHS [78] was used as base to develop this optimization algorithm. In the improvisation step, a roulette selection scheme was introduced where the selection probability of an individual is proportional to its fitness. This modification allows to consider the best individual in the evolutive process with a higher possibility. The independent minimization of the heat transfer area and the total pressure drop was performed where seven decision variables were considered: number of hot side layers, fin height, thickness, frequency and strip length, and hot and cold flow lengths. Discrete variables such as the number of hot side layers were handled as continuous, and their values were rounded up to evaluate the objective function. An adaptative penalization method was employed to transform the optimization problem from constrained to unconstrainted one. Algorithm parameters were as follows: bw_i minimum and maximum values were defined in terms of the optimization variables bounds, PAR_min = 0.1, PAR_max = 0.99, HMCR = 0.9, and HMS = 5. Results showed that the proposed method was effective for PFHEs design and also highlighted the benefits of using a self-adaptative penalization mechanism to resolve this design optimization problem.

Turgut et al. [229] studied the thermal design of STHE using a proposed enhanced version of the intelligent tuned harmony search algorithm (I-ITHS). In this method, a perturbation mechanism adopted from ABC is incorporated to enhance ITHS search capabilities. To improve the solution accuracy, OBL was used. Also, the chaos theory (i.e., Henon equation) was employed to increase the solution diversity inside the improvisation step. The total cost of the heat exchanger was minimized considering the capital investment and annual and total discounted operation costs. A detailed description of the STHE mathematical model was provided. The shell diameter, baffle spacing, number of passes, and outside tube diameter were the design variables. I-ITHS was successful to reduce the total cost employing less computational resources in comparison with previous studies. The authors concluded that I-ITHS showed a promising performance in engineering design problems, thus suggesting its use in the thermal system design and heat exchanger optimization.

3.3. Distillation Design and Optimization

Distillation is the main separation process in chemical engineering industry. A significant amount of fluids in the chemical industry (∼95%) is separated by this process, which is related to around 3% of global energy consumption [71]. This separation system has been studied for decades to select its optimal configuration and operating conditions to reduce costs, energy consumption, and environmental impacts [293–295]. The optimal design of distillation sequences for the separation of multi-component system is still one of the challenging design problems in process systems engineering [296]. The global and multi-objective optimization of a distillation system is known to be a problem involving discrete and continuous variables, as well as equality and inequality constraints that must be satisfied [297]. HS-based methods have been also applied in the distillation design.

Cabrera-Ruiz et al. [221] studied the performance of HS, SA, and GA in the design of three distillation sequences: conventional, dividing wall column, and thermally coupled reactive system with side stripper. It has been reported that the last two sequences can be more efficient (with up to 30% energy cost reduction) in comparison with the conventional scheme [297, 297]. The heat load was minimized considering both continuous and discrete variables. For the conventional distillation applied to a hydrodesulfurization process, the objective function was defined in terms of the feed temperature, reflux ratio, column pressure, feed stage, and total number of stages. For the dividing wall distillation column, the objective function was minimized with respect to the reflux ratio, column pressure, distillate flux, liquid and vapor interconnection flows, side fluxes, side stream tray location, first and last prefractioner tray location, feed stage, and total number of stages. This sequence was used for the purification of a mixture of alcohols. Finally, the thermally coupled reactive distillation sequence was used to obtain biodiesel from lauric acid and methanol. In this case, the optimization variables were the reflux ratio, column pressure, value and location of the interconnection flow, number of stages of the main column and stripper, first and last reaction tray location, stream vapor tray location, distillate and side fluxes, feed stage, and total number of stages. The optimization methods were coupled with the Aspen Plus simulator to carry out the minimization of objective functions of these processes. The bounds of the decision variables were defined based on preliminary calculations, and the HS parameters were PAR = 0.4, HMCR = 0.8, and HMS = 10. SA showed better results for the conventional and thermally coupled reactive distillation sequences in comparison to HS and GA, while HS was the best method for the design of dividing wall distillation column design.

The optimal design of heat integrated distillation systems is a relevant subject focused on thermal efficiency improvement via several techniques, which seek to effectively use heat from distillation columns to reduce the energy input requirements [71]. The synthesis of heat integrated systems requires optimizing three main interrelated aspects: the number and sequence of separation units, the individual design of the columns, and the heat exchanged network design [235] (see Figure 7 for an illustrative example). Several techniques have been used to solve this design problem including stochastic global optimization algorithms. Lashkajani et al. [235] developed a methodology to optimize the total annual cost of different heat integration distillation sequences involved in the olefin production using IHS [78]. First, the best separation sequence was found considering all flow rates, compositions, and several specifications for the whole separation network. This sequence was used as reference to evaluate the effect of heat integration in the total annual cost. The decision variables of the problem were the condenser and reboiler heat flows and temperatures, pressures, diameters, reflux ratios, cooling rates, actual plates, and heating rates. Three examples were analyzed to evaluate and compare the proposed methodology versus the results obtained by Isla and Cerda [301]. For all tested problems, IHS obtained the best performance. Additionally, an olefin production problem from a seven-hydrocarbon mixture was also solved where IHS results were compared with those obtained with GA. The final heat integrated sequence design proved to be an economical and competitive solution. The addition of heat integration techniques reduced the total annual cost by 44.64 and 35.37% using IHS and GA, respectively.

Figure 7 

Flowsheet of a heat integrated reaction distillation sequence corresponding to a minimum energy consumption configuration [300].

Boto et al. [239] applied HS and GA to minimize the energy and water consumption in a thermal cracking process. The mass flow rates of naphtha into the furnaces, steam dilution ratio, temperature, induction to turbines, and the distillate to feed ratio in the demethanizer column were defined as the optimization variables. These variables were selected using a sensitivity analysis to identify those with the highest impact on the energy and water reduction. The cracking process was constrained by the conditions of equipment performance, resource specifications, production requirements, and quality properties, which were considered in the optimization problem formulation. HS parameters used in the resolution of this problem design were HMCR = 0.7 and PAR = 0.5, and bw_i was defined in terms of a constant parameter and the lower and upper bounds of decision variables. HS achieved the best performance where 40% reduction of water and energy consumption was reported.

Sukpancharoen et al. [241] proposed a technique to synthesize a heat integrated distillation sequence using HS. The objective function was the total annual cost in terms of energy and the configuration of the separation sequence. A variance analysis was carried out to stablish the optimal parameters of the algorithm (i.e., PAR, HMCR, HMS, and iter_max). HS proved to be successful of solving large-scale mixed integer linear programming problems without any additional mathematical manipulation to suppress the non-convexities or to divide the problem into a set of subproblems.

Mansourzadeh et al. [244] implemented HS in a new multi-component isotope separation cascade code (MISCC-HS) for the calculation of optimal parameters of a centrifuge separation cascade for multi-component systems. Cascade optimization is considered a constrained problem where the optimal design is a configuration with a minimum total flow (or minimum number of gas centrifuges) for a given isotope concentration or a maximum separation capacity for a specific feed flow rate [244, 302]. The objective function was defined with respect to the stage cuts and feed location according to a defined minimum concentration in the distillate and a maximum concentration in the bottom. Constraint violations (i.e., out of bound concentrations) were managed using a penalty value. The proposed method was tested and compared with other strategies for the separation of uranium and xenon isotopes. The HS method proved to be suitable in the cascade separation design for mixtures.

Gao et al. [242] introduced a hybrid algorithm called cultural harmony search (CHS) to solve a real-time diesel mixture optimization problem, which is a complex high-dimensional problem with non-linear and linear constraints. The objective was to maximize the sum of flow rates of final products of this process subject to flow rate conservation and quality property constraints. The optimization variables were the distillates from each column, flow rates of the temporal blending tank, catalytic cracking, and cooking units. The variable bounds were calculated according to the initial feed flow rate of each column and the yield percentage. HS was implemented with the following parameters: PAR = 0.75, HMCR = 0.95, and HMS = 50. The performances of HS, CHS, improved CHS (ICHS), and simplex improved CHS (SICHS) were compared. The implementation of the simplex method to accelerate the convergence by a domain reduction in the decision variables was very effective to find better solutions to this optimization problem.

3.4. Life Cycle Analysis

Life cycle analysis (LCA) is a strategy to perform a comprehensive assessment of a product from raw material extraction to its production process, use, recycling, and disposal. LCA provides useful information for decision making in process design considering the environmental impacts, economic factors, energy consumption, and other relevant metrics [303]. For example, the construction of residential buildings can involve the optimization of several design objectives such as cost reduction, energy consumption, and environmental impact (e.g., solid wastes) minimization. In this regard, Fesanghary et al. [197] presented a multi-objective HS-based method to minimize the life cycle cost (LCC) and equivalent carbon dioxide (CO₂-eq) emissions in the construction of a residential building. The optimization variables were treated as discrete values, and only the market available concepts were considered in the optimization procedure (e.g., external walls, roofs, garage doors, windows, and floors). LCC analysis considered three phases: pre-use, use, and end of life of the building (see Figure 8). The first phase involved the extraction and processing of raw materials, component manufacturing, transport, and structure edification. In the use phase, CO₂ emissions related to energy consumption for lighting, heating, and cooling, as well as those related to maintenance were considered over a period of 25 years. Finally, the final phase comprised demolition of the building and transportation of residues to the landfill or recycle centers. Initial solution vectors were randomly generated in HS and then sent to the EnergyPlus simulator to calculate the environmental impact, energy consumption, and life cycle cost, which were used to evaluate the objective function in HS (see Figure 9). Results showed that it was only possible to reduce CO₂ emissions with an increase in LCC. Asadi [230] also minimized the LCC and CO₂-eq emissions for household design considering the whole structure life cycle (i.e., pre-use, use, and end of life phases). In this case, the transportation was not considered. Discrete decision variables were divided in structure-related parameters and climate regions. Results showed that there was a significant effect of climate regions in the final designs. Further research was suggested on the effect of different HVAC (heating, ventilation, and air conditioning) systems on optimal designs.

Figure 8 

Building life cycle analysis scheme [304].

Figure 9 

Illustration of schematic structure of the life cycle analysis [197].

Maleki et al. [236] used DHS [96] to minimize the total life cycle cost of a hybrid system composed of photovoltaic panels, a diesel generator, and a storage battery in a real-world application for a domestic house in an Iranian village. The total life cycle cost included the maintenance, capital, and total annual consumption costs of the diesel generator. This combinatorial optimization problem involved one continuous and two discrete decision variables: the diesel fuel consumption and the number of batteries and panels, respectively. The constraints on total load demand, charge/discharge battery power, diesel generator power, and solar power were considered. The parameters used for DHS were iter_max = 1,000, bw_min = 0.01, bw_max = 1, PAR_min = 0.1, PAR_max = 0.7, HMCR = 0.95, and HMS = 20. Six biodiesels were compared for the total life cycle cost analysis of the hybrid system. The best results were associated with the biodiesel produced from Norouzak seed oil. Other biofuels resulted in 9 to 43% increment of the total life cycle cost. The optimal configuration comprised 101 photovoltaic panels, 31 batteries, and 1 diesel generator. This study concluded that the proposed hybrid system was suitable and profitable to be implemented on that geographic location.

3.5. Design of Water Distribution Systems

The design of a water distribution network (WDN) involves the determination of the optimal size, location, and capacity of all operations involved in this process. WDN is composed of interconnected hydraulic components such as pumps, reservoirs, tanks, distribution pipes, and water treatment plants. The WDN design aims to establish the water supply efficiency from the source to its final destination. Several objectives such as optimal functioning assessment in terms of electricity cost and hydraulic behavior have been studied in WDN design. Note that the objective is to supply the correct amount of water with an adequate pressure and quality to the end user. These calculations are computationally expensive and time-consuming because of their non-linear models and high-dimension characteristics [305]. Generally, the liquid resource is pumped to higher elevations to meet specified delivery pressures to ensure an adequate water supply service. There are some factors to consider in order to achieve a constant water supply, for example, the water consumption is not uniform throughout the day, and water pumps cannot be activated instantly. In this regard, the water is usually stored in highly elevated tanks to ensure the supply. Therefore, the pumping operations imply significant costs, so the design of optimal energy-saving water supply schedules is a relevant subject [237].

Geem [216] introduced a cost optimization method for pumping in WDNs using HS. The optimization variables were the pipe diameter and pump size. HS was used as optimizer with the following parameters for this application: PAR ∈ (0.05 – 0.2), HMCR ∈ (0.95 – 0.99), and HMS ∈ (30 – 100). The maximum number of iterations was defined based on similar works, and a penalty function was used to handle the minimum pressure constraint. The objective function was defined in terms of three major cost-related terms: capital cost of pipes, capital cost of pumps, and energy cost for pumping. WDN example presented by Costa et al. [306] was solved using HS where a complete analysis of this network implicates 10¹² possible designs. Similar results were found by HS and those reported from Costa et al. [306] in terms of the optimum value. However, HS required less numerical effort over multiple trials to resolve this design problem.

Yoo et al. [232] presented a technique to estimate the optimal pipe diameter of an agricultural looped water irrigation network. In looped systems, the flow passing through each pipe is unknown, increasing the resolution complexity. The design cost was minimized according to the construction, pipe materials, and maintenance costs. Note that these terms are in the function of the pipe diameter. The optimization problem was solved as an unconstrained problem using the penalty method to handle the hydraulic restrictions. WDN design involved 356 pipes (solution variables) and 18 different diameters that were used in the cost evaluation. HS parameters were set as PAR = 0.01, HMCR = 0.97, and HMS = 30. Results showed that HS-based optimization decreased the costs by 9% in comparison with economic pipe calculations methods.

De Paola et al. [237] employed a HS multi-objective optimization (HSMO) approach to solve the pump scheduling problem. This stochastic technique was coupled with the known hydraulic solver EPANET 2.0 to evaluate the problem constraints and to calculate the performance of the proposed pumping schedules. Two objective functions were minimized: the energy cost and number of pump switches (i.e., turning on/off a previously non-operating/operating pump). The main restriction of the problem was imposed by the tank’s water level, which was bounded between storage head values in function of pump station flows and previous water levels. Additional constraints were also defined: deliverable pump flows and initial water level. Usually, they are the most common constraints in the pump scheduling problem. Since decision variables values are 0/1, the original pitch adjustment step was omitted, and an alternative strategy was used to generate the corresponding integer values.

Jung et al. [238] developed a hybrid algorithm based on HS and a pattern recognition method (HyHS) to perform WDN design. After 100 iterations, the solution vectors stored in HM were used as a training set to identify promising patters to find better solutions in future iterations. This algorithm was tested on the New York tunnel problem to minimize the cost with respect to the pipe sizes considering a minimum pressure constraint. Calculations showed that the integration of the pattern recognition method in the HS algorithm improved its robustness compared to IHS [78], IGA, and IGA + pattern recognition. The resolution of a more complex network using HyHS was suggested to evaluate its effectiveness under other challenging scenarios.

Yoo et al. [240] employed the HS algorithm to evaluate the optimal re-chlorination location and dosing in two real WDNs in Korea considering current and future conditions. In WDN design, the compliance with standard residual chloride concentration is necessary to avoid the appearance of diseases. The first study analyzed eight reservoirs that distributed 54,020 m³ of water through 85 and 773 km of transport and distribution pipes, respectively. For the second case, WDN with 3.8 km of pipes on 44 m³ surface was studied. The objective function for each problem was the mass minimization of additional chlorine that would be injected into the supply network nodes. The decision variables were the mass injection rates for all nodes. The problem was subject to the constraints imposed by the concentration of residual chlorine at each node, number of re-chlorination facilities, and conservation of mass and energy at each node. A penalty approach was employed to handle these constraints. A sensitivity analysis was carried out to determine the best values for the HS parameters, which were set as follows: bw = 0.01, PAR = 0.1, HMCR = 0.7, and HMS = 30. A maximum number of iterations was defined as the algorithm termination criterion. For both cases, the re-chlorination injection at four nodes in the network was the best solution. A comparison was made between HS and GA, and the former algorithm showed the best performance. Results showed that the WDN optimization complies with the residual chlorine concentration standards, so the proposed model proved to be useful for the WDN design.

3.6. Other Applications in Chemical Engineering

HS-based methods have been used to resolve a wide diversity of other global and multi-objective optimization problems related to process systems engineering. For example, Vasebi et al. [214] used HS to solve the combined heat and power economic dispatch problem. The total power and heat production cost was minimized in terms of heat and power production per unit subject to several constraints such as specific heat and power requirements. HS was applied using the following parameters: PAR = 0.5, HMCR = 0.85, and HMS = 6. A novel test system was presented and suggested as future reference to assess the algorithm performance in combined heat and power economic dispatch problems.

Zou et al. [219] also used the NGHS method to solve 39 chemical equation balance problems. The reactants and products’ stochiometric coefficients were the optimization variables to minimize a normalized constraint value related to the mass conservation equations. Real values were employed in the improvisation step and then converted to discrete stochiometric coefficient values. It was found that 3 of 39 chemical equations were challenging to solve. NGHS was reported as an efficient and robust method for chemical equation balancing.

Jiang et al. [227] proposed a HS variant called almost-parameter-free harmony search (APF-HS) to solve a groundwater pollution problem. The main objective of this problem was to identify the relevant characteristics of the pollution source(s) such as release events or fluxes. A weighted least-squares type objective function defined in terms of calculated and observed pollutant concentrations was used. One of the main features of the APF-HS method relies on the HMCR and PAR parameters which were dynamically adjusted according to algorithm iterations. Specific values of HMCR and PAR were used for each optimization variable. Regarding the bw parameter, its value was defined according to the minimum and maximum values in algorithm’s memory. A penalized objective function was used to handle the constraints. Comparative results between the proposed technique and IHS were reported.

Ma et al. [228] implemented the HS method in the Alopex-based evolutionary algorithm (AEA) [307] in order to improve its search mechanism. The proposed approach (HSAEA) was tested in two reaction kinetic parameter estimation problems. The first problem was related to the heavy oil thermal cracking three-lump model where eight parameters were estimated via the minimization of the average self-checking relative error. In the second case of study, ten parameters (i.e., the activation energies of the reaction kinetic equations and the pre-exponential factors) related to the reaction mechanism of mercury oxidation were determined. The average squared difference between the experimental and calculated mercury conversion was minimized. HSAEA parameters were bw_min = 1/(30 (UB_i − LB_i)), bw_max = 1/(20 (UB_i − LB_i)), PAR_min = 0.01, PAR_max = 0.99, and HMCR = 0.95. For both cases, the proposed HSAEA algorithm showed competitive results.

Dynamic optimization problems were solved by Fan et al. [308] using HS hybridized with differential evolution, which has co-evolutional control parameters. Dynamic optimization implies the resolution of problems described by differential equations that represent the change of design (target) variables with respect to time [309]. In particular, this hybrid HS method was assessed in the feeding-rate optimization of foreign protein production and an optimal control problem, thus concluding that this method outperformed the results reported by other authors. These results suggested that the application of HS-based method can be extended to the reactor network synthesis, parameter estimation of dynamic models, and optimal control [309].

Enikeeva et al. [251] used HS and gravitational search to solve a relevant problem in chemical rector design: the inverse problem of chemical kinetics. This problem is characterized for its wide search space, the presence of non-linear and non-differentiable functions, and also for being multi-modal. The optimization variables were a set of activation energies and pre-exponential coefficients for each reaction, as well as some additional parameters included in the reaction rate expressions. The inverse chemical kinetics problem can become a large-scale optimization problem due to complex chemical reactions involved in industrial operations. The proposed HS algorithm was evaluated using 20 benchmark functions including uni-modal and multi-modal objectives. The performance of this algorithm was assessed using the mean value and standard deviation of the decision variables. The former expressed the quality of the obtained solution, and the second one was associated to the algorithm stability. Two chemical operations were studied: the pre-reforming of propane into methane gas using a catalyst and the catalytic isomerization of the pentane-hexane fraction. HS was implemented with the following parameters: bw = 0.25, PAR = 0.95, HMCR = 0.95, and HMS = 100. These authors concluded that HS performed better in problems with low dimensionality.

Yadav et al. [250] studied the solution of a dynamic model for multiple effect evaporators using Fourier series and two metaheuristics, PSO and HS. Particularly, HS was used to optimize the Fourier series coefficients employed to solve the simultaneous non-linear ordinary differential equations of this problem. Fourteen variables were optimized with HS using the following algorithm parameters: iter_max = 1000, PAR = 0.1, HMCR = 0.9, and HMS = 30. The penalty method was used to handle the constraints for the liquor concentration and vapor temperature variables. Overall, the results showed that HS outperformed PSO.

4. Remarks on Challenges to Improve the Performance of HS for Chemical Engineering Optimization

HS is a metaheuristic algorithm that has drawn the attention of the scientific community since its initial proposal. In this review, the application of this algorithm in the chemical engineering field has been analyzed and discussed. First, a brief introduction covering general optimization concepts, classification, and characteristics of optimization techniques, previous reviews on HS, and motivating remarks for its use on chemical engineering were presented. Next, the fundaments of HS have been described in detail, along with popular variants reported in the literature. Recent developments on the algorithm, a wide variety of proposed variants, and multi-objective adaptations were also discussed. The applications of HS in chemical engineering were introduced to the readers in several sections covering single and multi-objective optimization problems with non-convex and non-linear behaviors subject to different types of constraints. This overview of the application of HS in chemical engineering indicated that this stochastic optimizer has been applied successfully in a wide diversity of problems. However, there are challenges and further directions of research that could improve the HS metaheuristic.

A great number of HS variants have been published over the years to enhance the method performance based on the modification/adaptation of its control parameters, structure, or via the hybridization with other metaheuristic operators. However, a limited number of these algorithm variants have been extensively assessed in the literature to characterize its reliability and efficiency to solve optimization problems with different dimensions, optimization variables, and constraints. A numerical characterization of existing and novel HS algorithms using a general testing framework could provide a reliable picture of their capabilities and limitations.

On the other hand, the resolution of constrained optimization problems with HS algorithms is commonly based on the penalty-based approach to handle the corresponding restrictions. This technique is strongly dependent on the penalty factor values, and further insights on the use of alternative constraint handling methods should be reported. Some studies have also highlighted the impact that the termination criterion could have on the algorithm robustness. Consequently, a generalized analysis of different stopping conditions for both global and multi-objective optimization problems should be performed. Until now, this analysis has been carried out for the application of HS in the parameter estimation of some thermodynamic models, Gibbs free energy minimization, and phase stability. The research on parameter-setting-free exploration and exploitation mechanisms could be promising to enhance HS capabilities to solve diverse chemical engineering optimization problems.

HS-based optimization methods have proved to be competitive numerical tools to resolve challenging optimization problems from chemical engineering. This stochastic method has a high potential to be applied in a wide window of other design problems like process control, supply chain design, and other emerging areas related to chemical engineering.

Data Availability

The data used to support the findings of this study are available from the corresponding authors upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

R. M. Fonseca-Pérez, O. Del-Mazo-Alvarado, and A. Meza-de-Luna were responsible for original draft preparation, conceptualization, methodology, and investigation. A. Bonilla-Petriciolet and Z. W. Geem were responsible for review and editing, investigation, supervision, and project administration.

Acknowledgments

Instituto Tecnológico de Aguascalientes provided financial support for this study.