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

Authors: R. M. Fonseca-Pérez; O. Del-Mazo-Alvarado; A. Meza-de-Luna; A. Bonilla-Petriciolet; Z. W. Geem
Journal: International Journal of Chemical Engineering
 (2022)

---

## 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.


---

## Body

## 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. Section2 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 as(1)Optimizefx,subject to(2)LBi≤xi≤UBi,i=1,2,…,n,where fx is the objective function, x=x1,x2,…,xn is the solution vector, xi is a decision variable subject to lower and upper bounds LBi and UBi, respectively, and n 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:(3)gix≤0,i=1,…,p,(4)hix=0,i=1,…,q,where q and p are the number of equality hix and inequality gix 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 values(5)HM=x11x21⋯xn1x12x22⋯xn2⋮⋮⋱⋮x1HMSx2HMS⋯xnHMS|fx1fx2⋮fxHMS.
Each candidate solution vector is evaluated and sorted in HM according to its value of objective function:fx1<fx2<⋯<fxHMS.Step 2.
Improvise a new harmony. In this stage, a new harmony (i.e., new solution vector)xnew=x1new,x2new,…,xnnew 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 xinew are obtained according to(6)xinew=xij∈xi1,xi2,…,xiHMS,with probability HMCR,randLBi,UBi,with probability1−HMCR,where 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 xinew is obtained via the memory consideration, the pitch adjustment procedure is performed according to the following operator:(7)xinew=xinew±rand0∗bw,with probability PAR∗HMCR,xinew,with probability1−PAR∗HMCR,where rand0 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 usingxinew; if its value is better than the worst stored in HM, fxnew<fxHMS, the new solution vector will substitute its place in HM.Step 4.
Iterative cycle. Steps2 and 3 are repeated until the specified termination criterion is satisfied.
Illustratively, Figure2 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 via(8)PARiter=PARmin+PARmax−PARminitermax∗iter,(9)bwiter=bwmaxexplnbwmin/bwmaxitermax∗iter,where 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 (itermax), 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 xinew 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.
Step 2: HS algorithmStep 2: GHS algorithm(1)for each i∈1,ndo(2) ifrand1≤HMCRthen(3) xinew=xij where j∈1,HMS(4) ifrand2≤PARthen(5) xinew=xinew±rand∗bw(6) end if(7) else(8) xinew=randLBi,UBi(9) end if(10) end forfor each i∈1,ndoifrand1≤HMCRthenxinew=xij where j∈1,HMSifrand2≤PARthenxinew=xkbest where k∈1,nend ifelsexinew=randLBi,UBiend ifend forSuperscript best corresponds to the index of the best harmony in HM.


#### 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 pm to avoid local solutions caused by premature convergence. Also, xnew replaces xHMS 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.
Step 2: NGHS algorithm(1)for each i∈1,ndo(2)  xR=2∗xibest−xiworst(3)  ifxR>UBithen(4)   xR=UBi(5)  else ifxR<LBithen(6)   xR=LBi(7)  end if(8)  xinew=xiworst+rand∗xR−xiworst(9)  ifrand≤pmthen(10)   xinew=randLBi,UBi(11)  end if(12) end forxR is the value of the decision variable with adaptative step. Superscript worst corresponds to the index of the worst harmony in HM.


#### 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:(10)bwiter=bwmax−bwmax−bwminitermax∗2⋅iter,if iter<itermax2,bwmin,if iter≥itermax2.


#### 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:(11)xinew=xij∈xi1,xi2,…,xiHMS+G0,1∗bw,with probability HMCR,randLBi,UBi,with probability1−HMCR,where a Gaussian distribution with mean 0 and standard deviation of 1 is used to calculate the random number G0,1. On the other hand, the pitch adjustment mechanism is given by(12)xinew=xibest+ϕ∗bw,with probability PAR∗HMCR,xinew,with probability1−PAR∗HMCR,where ϕ∈ [−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 xinew takes the value of xibest as in GHS. HMCR and PAR are dynamically adapted according to(13)HMCRiter=HMCRmin+HMCRmax−HMCRmin∗iteritermax,PARiter=PARmin+PARmax−PARmin∗iteritermax.


#### 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 by(14)HMCRiter∈randn1+n,11+n,PARiter=1−iter−1itermax,bwi=xih−xim,ifxjh≠xji,UBi−LBi∗exp−UBi−LBi∗iteritermax,otherwise,where superscripts h and m refer to random integers within [1, HMS]. The pitch adjustment mechanism of ESGHS is given by(15)xinew=xih+ϕ∗xim−xih.Ifxim=xih at any point, the following equation is used for the pitch adjustment:(16)xinew=xih+ϕ∗UBi−LBi∗exp−UBi−LBi∗iteritermax.In the random selection (1 − HMCR), the new harmony element is obtained as follows:(17)xinew=randGμi,σ,where a Gaussian distribution with mean μi=xibest and standard deviation σ=1−iter−1/itermax2 is used to calculate the random number randG.


#### 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 APiter is calculated at each iteration to decide if the new solution is accepted or rejected following the next criterion:(18)APiter=1,iffxnew≤fxHMS,fxHMS−fxbestfxnew−fxbest,iffxnew>fxHMS.


#### 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 Tables3–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.
AlgorithmRemarksReferenceHS∗Based on a detailed mathematical analysis of the HM variance, the bw parameter was defined as the population standard deviation.Mukhopadhyay et al. [89]HS∗Self-adaptive HS. The bw parameter is replaced into the new harmony calculations considering the minimum and maximum values in HM. Furthermore, low-discrepancy sequences are used to initialize HM in place of random number generation.Wang and Huang [90]HS∗Adaptive HS. Dynamical adjustment of PAR and HMCR by a probabilistic selection based on mean values of the HM matrix.Hasançebi et al. [91]PSF-HSParameter-setting-free HS. Three steps were proposed to avoid the tuning of HS parameters: (1) a random tuning where HM is filled with random solution vectors, (2) a rehearsal step where a certain amount of solution vectors is generated and an additional HMS size memory keeps track of harmony source, and (3) a performance step where HMCR and PAR are used for each decision variable.Geem and Sim [92]PHSProposed HS. The bw parameter is calculated using derivatives (finite difference method) and lower and upper variable bounds. The probability of pitch adjustment is modified. An improvement to PHS was also presented where two bw definitions are used in the improvisation step. New decision probabilities are employed.Jaberipour and Khorram [93]EHSExplorative HS. The standard deviation of the current population defines the value of bw. Penalty strategy is used for constraint handling.Das et al. [94]NGHSA NGHS adaption to solve large-scale knapsack problems. The removal of HMCR and PAR is proposed to introduce a genetic mutation concept. A strategy from PSO is used to guide the worst solutions to global-best surroundings.Zou et al. [95]DHSDiscrete HS. It is designed to handle discrete variables. PAR and bw are set considering the defined bounds as well as the current and maximum number of iterations.Askarzadeh [96]NSHSNovel self-adaptive HS. HMCR and bw are fixed according to the problem dimension and decision variable bounds.Luo [97]LAHSLearning automata-based HS. Smart parameter setting of HMCR, PAR, and bw through a learning automata implementation [98].Enayatifar et al. [99]DSHSDynamic self-adaptive HS. It is based on the implementation of the best-to-worst ratio (BtW) concept in HM. PAR is dynamically adjusted in terms of BtW. Two bandwidth parameters are used, active bw and dynamic bandwidth value dbw. These parameters are calculated for each decision variable and are dynamically adjusted according to its standard deviation in all the HM.Kattan and Abdullah [100]ABHSAdaptative binary HS. New harmonies are created based on single or multiple vectors in HM. PAR is modified according to the best HM vector. HMCR and PAR are obtained using several mechanisms: linear decrease, linear and non-linear increase, and random, where all are calculated as a function of current and total iterations. ABHS can generate one or more harmonies at each generation. HM update is performed in serial or parallel according to the number of generated new candidates.Wang et al. [101]ABHSAdjustable bandwidth HS. Minimum and maximum values for PAR and HMCR are defined. Additionally, a saturation parameter (SP) is proposed to control bw parameter. An alternative version of ABHS with changes in bw calculation is also presented.Contreras et al. [102]GDHSGlobal dynamic HS. It is based on a dynamic change in bw, PAR, and HMCR parameters, as well as in the domain of HM variables.Khalili et al. [103]PAHSParameter adaptive HS. Four dynamic change strategies for PAR and HMCR in terms of the number of iterations were implemented.Kumar et al. [104]SFHSSelf-regulated fretwidth HS. It considers changes in bw (also called fretwidth parameter) at each iteration based on several criteria.Amaya et al. [105]SBHSSimplified binary HS. It focuses on changes to the improvisation step. The method is parameter free (i.e., PAR and bw were removed), and HMCR is dynamically adjusted. It incorporates greedy local search. SBHS was used to solve large-scale knapsack problems.Kong et al. [106]HS∗Dynamic adjustment of bw. For a particular decision variable, bw is defined according to the current iteration and its maximum and minimum values.Kalivarapu et al. [107]GGHSGaussian GHS. HM is initialized based on a Gaussian distribution and it is adjusted considering a dynamic bw. The best harmony is also adjusted using a dynamic standard deviation.Keshtegar and Sadeq [108]ImHSImproved HS. The bw parameter decrements according to an exponential function in terms of the number of iterations. A selection scheme based on a proposed parameter called rate of improvement adjustment is applied.Portilla-Flores et al. [109]SSaHSSemi-self-adaptive HS. The bw parameter is updated according to a self-adaptative rule, which considers the information about the current population. Two new parameters are stablished by an empirical analysis.Zhao et al. [110]PSFHSModified parameter-setting-free HS. The number of harmonies in HM is modified by a positive integer coefficient. PAR and HMCR domains are redefined considering the minimum and maximum values, respectively. Recommended limits for HMCR and PAR were provided.Shaqfa and Orbán [111]BWM_HSBest-worst-mean HS. It considers the replacement of HMCR for three novel criteria to enhance exploration/exploitation balance. A self-adaptative probabilistic factor is incorporated to choose between rules.Talaei et al. [112]FHSFuzzy HS. A fuzzy system is used to handle HMCR and PAR throughout iterations.Peraza et al. [113]EHSEnhanced HS. It considers dynamic adjustment of PAR and bw. PSO strategy is implemented in new harmony generation.Niu et al. [114]AHSAdaptative HS. Parameters HMCR and bw are defined as function of current and total number of iterations.Yücel et al. [115]SIHSSelf-adaptative improved HS. It considers the dynamic modification of bw starting with a large value and adjusting it at each generation.Houari et al. [116]∗Variant of HS.Table 4
Harmony search variants based on the algorithm structure modification.
AlgorithmRemarksReferenceHS∗Algorithm for structural optimization. It implies the initialization of HM memory using the finite element method considering problem’s constraints. This approach is also used during the improvisation step to assess new harmonies.Lee and Geem [117]HS∗It was proposed to solve discrete variable optimization problems employing a stochastic derivative formulation based on several probability criteria using HS parameters. It was tested on benchmark and real-world problems.Geem [118]HS∗HS with SA concept. It is based on the modification of PAR and bw. Starting from a maximum initial value, the parameters are slowly reduced imitating the cooling concept of the SA algorithm.Taherinejad [119]IHSImproved HS. It considers the dynamic modification of PAR based on dispersed particle swarm optimization (DPSO) [120]. The information index, which relates the difference of the particles according to its fitness value, is incorporated.Coelho and Bernert [121]DHSDifferential HS. DE mutation strategy is implemented in the improvisation step.Chakraborty et al. [122]PSHSParticle swarm HS. A concept of particle swarm optimization is used in the improvisation step.Geem [123]CHSIt uses chaotic maps (e.g., logistic, tent, Gauss, and so on) to substitute generally employed random numbers in HS operators. The algorithm was tested in multi-modal benchmark functions.Alatas [124]GHS + LEMMachine learning tools are employed to create new populations in HS according to the learnable evolution model (LEM) [125]. Also, an evolutionary related method called the Darwinian is used to decide which members of the population are better to fulfil certain tasks.Cobos et al. [126]ITHSIntelligent tuned HS. HM is divided in two groups in order to balance exploration and exploitation. The first one contains those solutions with an objective value less or equal to the mean value in HM, and the second group includes the rest of solution vectors. PAR is dynamically adjusted.Yadav et al. [127]NDHSA tournament selection rule is used to choose the candidate harmony from the HM. The solution vectors with a better objective function value will be more likely to be used to generate new solutions. Also, the PAR and bw parameters are dynamically adjusted.Chen et al. [128]HSTLHS based on teaching-learning. It was proposed to solve complex high dimension optimization problems. It considers the introduction of a teaching-learning concept into the improvisation step and a dynamical adjustment of all parameters (i.e., HMCR, PAR, bw, and TLP) in terms of generations.Tuo et al. [129]IGHSImproved GHS. The opposition-based learning (OBL) [130] was used as an alternative initialization procedure to improve the quality of the initial memory vectors. In addition, a mutation strategy was implemented in place of the GHS pitch adjustment. The random selection step was substituted for a modification used in artificial bee colony (ABC) algorithm [131].Xiang et al. [132]IGHSIntelligent global HS. In the improvisation step, the new harmony mimics an element of the best solution vector stored in HM.Valian et al. [133]GSHSGeometric selective HS. Two harmonies are combined to create theith element of the new harmony using a tournament selection strategy. A mutation operator is implemented using PAR.Castelli et al. [134]OHSOpposition-based HS. OBL strategy is incorporated in the HM initialization and update steps.Banerjee et al. [135]HS∗Several mutation strategies are incorporated in the random consideration process (i.e., random, boundary, non-uniform, MPT, power, and polynomial).Hasan et al. [136]QOHSQuasi-oppositional HS. It considers the introduction of the quasi-oppositional-based learning strategy in HS and the application of the quasi-opposition concept for HM initialization. PAR and bw are obtained according to current generation. It also implies the improvisation of new harmonies and populations on quasi-oppositional basis and jumping probability according to search space bounds.Shiva et al. [137]iHSIsland HS. It considers the incorporation of the island model [138] in the HM memory.Al-Betar et al. [139]DH/bestRandom initialization of HM is replaced by a new technique, which reduces randomness. The conventional pitch adjustment method is substituted by a DE mutation strategy.Abedinpourshotorban et al. [140]MHSModified HS. The bw parameter is removed. An intersect operation considering PAR is included, and after each improvisation, the algorithm performs a cellular local search.Yi et al. [141]LHSImproved HS. It considers three main modifications: an adaptative global pitch adjustment, incorporation of an OBL technique, and a competition selection strategy.Ouyang et al. [142]SRHSSelective refining HS. The best harmony is considered as a representation of the current experience and knowledge of the HM. A tournament selection scheme is used. Instead of taking decisions based only on the value of the objective function, the value of every decision variable is considered.Shabani et al. [143]GOGHSGlobal HS with generalized OBL. Relevant information from the global-best harmony is extracted to create new candidate solutions. OBL is employed.Guo et al. [144]TPHSTwo-phase HS. The algorithm is divided in two phases, the first is highly focused in exploration by the implementation of the catastrophic mutation concept used in GA, and the second phase uses local search to improve exploitation.Assad and Deep [145]DSAHSEnhanced HS with dual strategies and adaptive parameters. An improved pitch adjustment operation with dual strategies is incorporated to enhance the exploitation ability of the algorithm. The operation may use a best harmony or a random harmony-guided search strategy according to a probability decision. Furthermore, the parameters of the algorithm are self-adjusted according to the new candidate fitness.Wang et al. [146]HS-SAInspired by the SA concept, a new parameter called temperature is used so that inferior harmonies could also be considered. Two additional SA parameters are employed.Assad and Deep [147]DHSDifferential-based HS. It considers the introduction of a position indicator to enable continuous HS to solve the discrete job shop scheduling problem. It also includes the incorporation of the DE mutation concept into the pitch adjustment step.Zhao et al. [148]CIHSChaotic improved HS. The concept of chaos is introduced into improved HS using logistic map to calculate PAR and bw. The best solutions are refined using a chaotic local search. New harmonies emulate the best harmony employing a PSO concept.Shafaati and Mojallali [149]ABHSAdaptative HS with best-based search strategy. It employs a best-based search strategy, which makes use of the global-best solution to enhance the search. HMCR and PAR are self-tuned based on information from the evolutionary search. The bw parameter is not considered in this algorithm.Guo et al. [150]DANGHSDynamic adjusting NGHS. Several dynamic parameter adjustment strategies to calculate the genetic mutation probability were used: linear, threshold, exponential, and concave-convex.Chiu et al. [151]AHSAmended HS. Two perturbation strategies are used: stochastic differential and OBL. The bw parameter is dynamically adjusted in terms of the mean value of a particular decision variable. A global dimension selection to balance exploration and exploitation is used.Ouyang et al. [152]MHS-PCLSParallel chaotic local search enhanced HS. It considers an initial point modification in PCLS to improve the performance of HS. It also uses MHS intersect mutation operation.Yi et al. [57]IHSGlobal learning and a modified random selection operation are proposed. A new improvement scheme based on the iteration modes (i.e., dimension-to-dimension, stochastic selection, vector, and matrix) is proposed. A total of four IHS variants were developed.Ouyang et al. [153]PHSβ-HC and imp. PHSβ-HCBest polynomial HS algorithm with bestβ-hill climbing. Two variants of the proposed method were introduced: PHSβ-HC and Imp. PHSβ-HC. The random selection step in HS is replaced by roulette wheel and tournament selection. A highly disruptive polynomial mutation replaces the pitch adjustment and random selection steps.Doush and Santos [154]EHS-CRPEnhanced HS with circular region perturbation. A global and local dimension selection is used by considering the global-best solution vector, and also local partial dimension variables are randomly selected. The selection operation is based on several PSO concepts. The circular region perturbation method is employed to improve search efficiency.Wu et al. [155]IMGHSAIntersect mutation global HS algorithm. Two methods were combined: intersect mutation and global HS. It considers the introduction of a new operator to reinforce exploration. HM is divided in two and a new operator considers both parts to generate new harmonies.Gholami et al. [156]SLHSSequence learning HS. It was adapted to solve the combinatorial flexible process planning problem. The optimal sequence is stablished by finding a proper successor for each operation.Luo [157]AHS-HCMAmeliorated HS with hybrid convergence mechanism. New harmonies are calculated through a convergence coefficient and concepts of PSO and DE based on the number of iterations. Non-linear dynamic domain in terms of current and maximum number of iterations when 1-HMCR.Zhu and Tang [158]MHSAModified HS algorithm. It considers the HM initialization by means of opposition numbers. Two new parameters inspired by the GWO method were introduced to calculate new harmonies. HMCR is calculated based on a Gaussian distribution. PAR is formulated as a non-linear exponential function.Gupta [159]AHS-DE-OBLAdaptative HS using DE and OBL. It includes the generation of new harmonies with OBL approach and the dynamic adjustment of HMCR, PAR, and search domain according to current and maximum number of iterations. The best and worst harmonies are used within a DE concept to calculate bw.Kang et al. [160]∗Variant of HS.Table 5
Harmony search variants based on hybridization with other metaheuristics.
AlgorithmRemarksReferenceSA-HSNelder–Mead simplex algorithm HS. The simplex method is incorporated in HM update using then best HM vectors + 1. The n best HM harmonies are kept in memory and the next vector is updated with the simplex solution. The HS improvisation stage considers these vectors along with the remaining harmonies in HM.Jang et al. [161]HHSAHybrid HS algorithm. Sequential quadratic programming (SQP) is used to enhance the local search performance. Two strategies are proposed: use of best harmonies in HM as initial points for the local search; performing a local search for each improvised solution vector. The strategy selection is based on a probabilistic parameter.Fesanghary et al. [162]HTHSAHybrid Taguchi HS algorithm. Two optimization phases are considered: (1) the Taguchi method is used to fine-tune decision variables intervals and to generate initial solutions and (2) HS is used to perform the global search.Yildiz [163]HS-DEDE is used to fine-tune HM elements considering all solution vectors as the individuals of the population.Gao et al. [164]HS-CSAThe clonal selection algorithm is used to improve the elements of the HM.Wang et al. [165]HS-solverIntegration of HS with non-linear excel solver algorithm. Two optimization approaches were proposed. The first one consists of the use of the local method once the HS solution is obtained. In the second approach, the algorithms operate simultaneously and new HS solutions are optimized with the local method at each iteration.Ayvaz et al. [166]HMGHSHybrid modified GHS. It was proposed to solve the block permutation flow shop scheduling problem. Continuous variables were coded to binary using the largest position value rule. HM initialization is based on problem’s heuristic rule. It includes a local search feature.Wang et al. [167]DMS-PSO-HSHS is incorporated into each subswarm of the dynamic multi-swarm particle swarm optimizer algorithm [168], in order to benefit from the particular exploration capabilities of the PSO technique and the exploitation abilities of HS.Zhao et al. [169]DCHSSADiscrete chaotic HS-based simulated annealing. New solutions of the SA algorithm are generated using HS parameters (HMCR and PAR). Chaotic mapping is used to generate new solutions when the probability is within HMCR.Askarzadeh [170]HS/BAHybridization of HS with the bat algorithm (BA) [171] where the mutation operator is modified by applying the concept of pitch adjustment to improve convergence.Wang and Guo [172]HS/FAHS and firefly algorithm (FA) [173]. Hybrid method that combines HS for exploration and FA for exploitation. It also incorporates HMCR and PAR to perform a pitch adjustment in FA mutation step.Guo et al. [174]AntHSHS and ACO. It considers the addition of two ACO concepts to HS: pheromone values and heuristics.Amini and Ghaderi [175]GHSACOGHS with ACO. A probabilistic factor is added to switch between using GHS or ACO for the improvisation step. The use of ACO was proposed to enhance exploitation.Fouad et al. [176]CPSO-MHSCooperative particle swarm optimizer-modified HS. It is a hybrid algorithm comprising a PSO variant (CPSO) and a modified HS version. MHS is employed for exploration and CPSO for exploitation. Note that PAR and bw parameters were removed using CPSO for local search.Zhang and Li [177]HS/CSHS with CS. It applies HS strategy in CS in order to increase diversity of cuckoo populations in the optimization process. Specifically, the HS improvisation strategy is used in the mutation stage of CS.Wang et al. [178]HSTLBOHS with teaching learning-based optimization [179]. A dynamic selection strategy is used to determine which algorithm will perform the search for a given iteration is used. It also includes a dynamic change of PAR and bw. In the improvisation step, a dynamic strategy considering elements of the worst harmony is implemented.Tuo et al. [180]MBHSMine blast HS. It employs MBA and HS for the exploration and exploitation phases, respectively. HS uses concepts from IHS [78] and PSF-HS to calculate PAR and HMCR in terms of iterations. A new iteration-dependent dynamic harmonic memory is used.Sadollah et al. [181]HS-GWOHS with grey wolf optimizer. GWO is used to control PAR and bw parameters during the search, and OBL enhances the improvisation step.Alomoush et al. [182]H-TSHS with TS. The neighborhood search capabilities of TS are exploited after HS improvisation step. It was applied to the dynamic parallel row ordering problem.Gong et al. [183]HHSDEHybrid harmony DE algorithm. It improved mutation operation in DE making use of HS improvisation step. HMS is dynamically updated according to current iteration. Use of a self-parameter adjustment strategy to control method’s parameters.Fu et al. [184]HABC-HSA-MAHybrid ABC and HS algorithm-based metaheuristic approach. These algorithms are used sequentially, and the final output of HS optimization process is the input for the ABC method.Jayalakshmi et al. [185]HS∗HS combined with L-BFGS. After generating HS initial population, a local search using the quasi-Newton method is performed for each harmony, keeping record of the best value in HM.Sorokovikov and Gornov [186]UKF-HSUnscented Kalman filter HS. It was developed to solve the finite-element model updating using single and multi-objective optimization.Naranjo-Pérez et al. [187]AHHSAdaptive-hybrid HS. The Java algorithm is used in HS improvisation step. HMCR and bw are randomly defined at the beginning of the algorithm.Yücel et al. [115]HSCAHSHybrid sine-cosine approach HS. HS is used in parallel with SCA to improve solution quality. HS generates a new population seeking to improve the worst harmonic in HM.Singh and Kaur [188]RDEGA-HSReward population-based differential genetic HS. It combines two algorithms proposed by the authors. The first, reward-population HS (RHS), works with a set of independent subharmonies and a reward harmony. The second, differential genetic HS (DEGA-HS), removes the bw parameter and employs two DE mutation strategies and a concept of the parthenogenetic algorithm in HS improvisation stage.Zhang et al. [189]∗Variant of HS.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 as(19)Optimizef1x,f2x,…,fmx,where 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.


## 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 as(1)Optimizefx,subject to(2)LBi≤xi≤UBi,i=1,2,…,n,where fx is the objective function, x=x1,x2,…,xn is the solution vector, xi is a decision variable subject to lower and upper bounds LBi and UBi, respectively, and n 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:(3)gix≤0,i=1,…,p,(4)hix=0,i=1,…,q,where q and p are the number of equality hix and inequality gix 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 values(5)HM=x11x21⋯xn1x12x22⋯xn2⋮⋮⋱⋮x1HMSx2HMS⋯xnHMS|fx1fx2⋮fxHMS.
Each candidate solution vector is evaluated and sorted in HM according to its value of objective function:fx1<fx2<⋯<fxHMS.Step 2.
Improvise a new harmony. In this stage, a new harmony (i.e., new solution vector)xnew=x1new,x2new,…,xnnew 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 xinew are obtained according to(6)xinew=xij∈xi1,xi2,…,xiHMS,with probability HMCR,randLBi,UBi,with probability1−HMCR,where 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 xinew is obtained via the memory consideration, the pitch adjustment procedure is performed according to the following operator:(7)xinew=xinew±rand0∗bw,with probability PAR∗HMCR,xinew,with probability1−PAR∗HMCR,where rand0 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 usingxinew; if its value is better than the worst stored in HM, fxnew<fxHMS, the new solution vector will substitute its place in HM.Step 4.
Iterative cycle. Steps2 and 3 are repeated until the specified termination criterion is satisfied.
Illustratively, Figure2 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 via(8)PARiter=PARmin+PARmax−PARminitermax∗iter,(9)bwiter=bwmaxexplnbwmin/bwmaxitermax∗iter,where 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 (itermax), 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 xinew 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.
Step 2: HS algorithmStep 2: GHS algorithm(1)for each i∈1,ndo(2) ifrand1≤HMCRthen(3) xinew=xij where j∈1,HMS(4) ifrand2≤PARthen(5) xinew=xinew±rand∗bw(6) end if(7) else(8) xinew=randLBi,UBi(9) end if(10) end forfor each i∈1,ndoifrand1≤HMCRthenxinew=xij where j∈1,HMSifrand2≤PARthenxinew=xkbest where k∈1,nend ifelsexinew=randLBi,UBiend ifend forSuperscript best corresponds to the index of the best harmony in HM.


### 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 pm to avoid local solutions caused by premature convergence. Also, xnew replaces xHMS 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.
Step 2: NGHS algorithm(1)for each i∈1,ndo(2)  xR=2∗xibest−xiworst(3)  ifxR>UBithen(4)   xR=UBi(5)  else ifxR<LBithen(6)   xR=LBi(7)  end if(8)  xinew=xiworst+rand∗xR−xiworst(9)  ifrand≤pmthen(10)   xinew=randLBi,UBi(11)  end if(12) end forxR is the value of the decision variable with adaptative step. Superscript worst corresponds to the index of the worst harmony in HM.


### 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:(10)bwiter=bwmax−bwmax−bwminitermax∗2⋅iter,if iter<itermax2,bwmin,if iter≥itermax2.


### 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:(11)xinew=xij∈xi1,xi2,…,xiHMS+G0,1∗bw,with probability HMCR,randLBi,UBi,with probability1−HMCR,where a Gaussian distribution with mean 0 and standard deviation of 1 is used to calculate the random number G0,1. On the other hand, the pitch adjustment mechanism is given by(12)xinew=xibest+ϕ∗bw,with probability PAR∗HMCR,xinew,with probability1−PAR∗HMCR,where ϕ∈ [−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 xinew takes the value of xibest as in GHS. HMCR and PAR are dynamically adapted according to(13)HMCRiter=HMCRmin+HMCRmax−HMCRmin∗iteritermax,PARiter=PARmin+PARmax−PARmin∗iteritermax.


### 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 by(14)HMCRiter∈randn1+n,11+n,PARiter=1−iter−1itermax,bwi=xih−xim,ifxjh≠xji,UBi−LBi∗exp−UBi−LBi∗iteritermax,otherwise,where superscripts h and m refer to random integers within [1, HMS]. The pitch adjustment mechanism of ESGHS is given by(15)xinew=xih+ϕ∗xim−xih.Ifxim=xih at any point, the following equation is used for the pitch adjustment:(16)xinew=xih+ϕ∗UBi−LBi∗exp−UBi−LBi∗iteritermax.In the random selection (1 − HMCR), the new harmony element is obtained as follows:(17)xinew=randGμi,σ,where a Gaussian distribution with mean μi=xibest and standard deviation σ=1−iter−1/itermax2 is used to calculate the random number randG.


### 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 APiter is calculated at each iteration to decide if the new solution is accepted or rejected following the next criterion:(18)APiter=1,iffxnew≤fxHMS,fxHMS−fxbestfxnew−fxbest,iffxnew>fxHMS.


### 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 Tables3–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.
AlgorithmRemarksReferenceHS∗Based on a detailed mathematical analysis of the HM variance, the bw parameter was defined as the population standard deviation.Mukhopadhyay et al. [89]HS∗Self-adaptive HS. The bw parameter is replaced into the new harmony calculations considering the minimum and maximum values in HM. Furthermore, low-discrepancy sequences are used to initialize HM in place of random number generation.Wang and Huang [90]HS∗Adaptive HS. Dynamical adjustment of PAR and HMCR by a probabilistic selection based on mean values of the HM matrix.Hasançebi et al. [91]PSF-HSParameter-setting-free HS. Three steps were proposed to avoid the tuning of HS parameters: (1) a random tuning where HM is filled with random solution vectors, (2) a rehearsal step where a certain amount of solution vectors is generated and an additional HMS size memory keeps track of harmony source, and (3) a performance step where HMCR and PAR are used for each decision variable.Geem and Sim [92]PHSProposed HS. The bw parameter is calculated using derivatives (finite difference method) and lower and upper variable bounds. The probability of pitch adjustment is modified. An improvement to PHS was also presented where two bw definitions are used in the improvisation step. New decision probabilities are employed.Jaberipour and Khorram [93]EHSExplorative HS. The standard deviation of the current population defines the value of bw. Penalty strategy is used for constraint handling.Das et al. [94]NGHSA NGHS adaption to solve large-scale knapsack problems. The removal of HMCR and PAR is proposed to introduce a genetic mutation concept. A strategy from PSO is used to guide the worst solutions to global-best surroundings.Zou et al. [95]DHSDiscrete HS. It is designed to handle discrete variables. PAR and bw are set considering the defined bounds as well as the current and maximum number of iterations.Askarzadeh [96]NSHSNovel self-adaptive HS. HMCR and bw are fixed according to the problem dimension and decision variable bounds.Luo [97]LAHSLearning automata-based HS. Smart parameter setting of HMCR, PAR, and bw through a learning automata implementation [98].Enayatifar et al. [99]DSHSDynamic self-adaptive HS. It is based on the implementation of the best-to-worst ratio (BtW) concept in HM. PAR is dynamically adjusted in terms of BtW. Two bandwidth parameters are used, active bw and dynamic bandwidth value dbw. These parameters are calculated for each decision variable and are dynamically adjusted according to its standard deviation in all the HM.Kattan and Abdullah [100]ABHSAdaptative binary HS. New harmonies are created based on single or multiple vectors in HM. PAR is modified according to the best HM vector. HMCR and PAR are obtained using several mechanisms: linear decrease, linear and non-linear increase, and random, where all are calculated as a function of current and total iterations. ABHS can generate one or more harmonies at each generation. HM update is performed in serial or parallel according to the number of generated new candidates.Wang et al. [101]ABHSAdjustable bandwidth HS. Minimum and maximum values for PAR and HMCR are defined. Additionally, a saturation parameter (SP) is proposed to control bw parameter. An alternative version of ABHS with changes in bw calculation is also presented.Contreras et al. [102]GDHSGlobal dynamic HS. It is based on a dynamic change in bw, PAR, and HMCR parameters, as well as in the domain of HM variables.Khalili et al. [103]PAHSParameter adaptive HS. Four dynamic change strategies for PAR and HMCR in terms of the number of iterations were implemented.Kumar et al. [104]SFHSSelf-regulated fretwidth HS. It considers changes in bw (also called fretwidth parameter) at each iteration based on several criteria.Amaya et al. [105]SBHSSimplified binary HS. It focuses on changes to the improvisation step. The method is parameter free (i.e., PAR and bw were removed), and HMCR is dynamically adjusted. It incorporates greedy local search. SBHS was used to solve large-scale knapsack problems.Kong et al. [106]HS∗Dynamic adjustment of bw. For a particular decision variable, bw is defined according to the current iteration and its maximum and minimum values.Kalivarapu et al. [107]GGHSGaussian GHS. HM is initialized based on a Gaussian distribution and it is adjusted considering a dynamic bw. The best harmony is also adjusted using a dynamic standard deviation.Keshtegar and Sadeq [108]ImHSImproved HS. The bw parameter decrements according to an exponential function in terms of the number of iterations. A selection scheme based on a proposed parameter called rate of improvement adjustment is applied.Portilla-Flores et al. [109]SSaHSSemi-self-adaptive HS. The bw parameter is updated according to a self-adaptative rule, which considers the information about the current population. Two new parameters are stablished by an empirical analysis.Zhao et al. [110]PSFHSModified parameter-setting-free HS. The number of harmonies in HM is modified by a positive integer coefficient. PAR and HMCR domains are redefined considering the minimum and maximum values, respectively. Recommended limits for HMCR and PAR were provided.Shaqfa and Orbán [111]BWM_HSBest-worst-mean HS. It considers the replacement of HMCR for three novel criteria to enhance exploration/exploitation balance. A self-adaptative probabilistic factor is incorporated to choose between rules.Talaei et al. [112]FHSFuzzy HS. A fuzzy system is used to handle HMCR and PAR throughout iterations.Peraza et al. [113]EHSEnhanced HS. It considers dynamic adjustment of PAR and bw. PSO strategy is implemented in new harmony generation.Niu et al. [114]AHSAdaptative HS. Parameters HMCR and bw are defined as function of current and total number of iterations.Yücel et al. [115]SIHSSelf-adaptative improved HS. It considers the dynamic modification of bw starting with a large value and adjusting it at each generation.Houari et al. [116]∗Variant of HS.Table 4
Harmony search variants based on the algorithm structure modification.
AlgorithmRemarksReferenceHS∗Algorithm for structural optimization. It implies the initialization of HM memory using the finite element method considering problem’s constraints. This approach is also used during the improvisation step to assess new harmonies.Lee and Geem [117]HS∗It was proposed to solve discrete variable optimization problems employing a stochastic derivative formulation based on several probability criteria using HS parameters. It was tested on benchmark and real-world problems.Geem [118]HS∗HS with SA concept. It is based on the modification of PAR and bw. Starting from a maximum initial value, the parameters are slowly reduced imitating the cooling concept of the SA algorithm.Taherinejad [119]IHSImproved HS. It considers the dynamic modification of PAR based on dispersed particle swarm optimization (DPSO) [120]. The information index, which relates the difference of the particles according to its fitness value, is incorporated.Coelho and Bernert [121]DHSDifferential HS. DE mutation strategy is implemented in the improvisation step.Chakraborty et al. [122]PSHSParticle swarm HS. A concept of particle swarm optimization is used in the improvisation step.Geem [123]CHSIt uses chaotic maps (e.g., logistic, tent, Gauss, and so on) to substitute generally employed random numbers in HS operators. The algorithm was tested in multi-modal benchmark functions.Alatas [124]GHS + LEMMachine learning tools are employed to create new populations in HS according to the learnable evolution model (LEM) [125]. Also, an evolutionary related method called the Darwinian is used to decide which members of the population are better to fulfil certain tasks.Cobos et al. [126]ITHSIntelligent tuned HS. HM is divided in two groups in order to balance exploration and exploitation. The first one contains those solutions with an objective value less or equal to the mean value in HM, and the second group includes the rest of solution vectors. PAR is dynamically adjusted.Yadav et al. [127]NDHSA tournament selection rule is used to choose the candidate harmony from the HM. The solution vectors with a better objective function value will be more likely to be used to generate new solutions. Also, the PAR and bw parameters are dynamically adjusted.Chen et al. [128]HSTLHS based on teaching-learning. It was proposed to solve complex high dimension optimization problems. It considers the introduction of a teaching-learning concept into the improvisation step and a dynamical adjustment of all parameters (i.e., HMCR, PAR, bw, and TLP) in terms of generations.Tuo et al. [129]IGHSImproved GHS. The opposition-based learning (OBL) [130] was used as an alternative initialization procedure to improve the quality of the initial memory vectors. In addition, a mutation strategy was implemented in place of the GHS pitch adjustment. The random selection step was substituted for a modification used in artificial bee colony (ABC) algorithm [131].Xiang et al. [132]IGHSIntelligent global HS. In the improvisation step, the new harmony mimics an element of the best solution vector stored in HM.Valian et al. [133]GSHSGeometric selective HS. Two harmonies are combined to create theith element of the new harmony using a tournament selection strategy. A mutation operator is implemented using PAR.Castelli et al. [134]OHSOpposition-based HS. OBL strategy is incorporated in the HM initialization and update steps.Banerjee et al. [135]HS∗Several mutation strategies are incorporated in the random consideration process (i.e., random, boundary, non-uniform, MPT, power, and polynomial).Hasan et al. [136]QOHSQuasi-oppositional HS. It considers the introduction of the quasi-oppositional-based learning strategy in HS and the application of the quasi-opposition concept for HM initialization. PAR and bw are obtained according to current generation. It also implies the improvisation of new harmonies and populations on quasi-oppositional basis and jumping probability according to search space bounds.Shiva et al. [137]iHSIsland HS. It considers the incorporation of the island model [138] in the HM memory.Al-Betar et al. [139]DH/bestRandom initialization of HM is replaced by a new technique, which reduces randomness. The conventional pitch adjustment method is substituted by a DE mutation strategy.Abedinpourshotorban et al. [140]MHSModified HS. The bw parameter is removed. An intersect operation considering PAR is included, and after each improvisation, the algorithm performs a cellular local search.Yi et al. [141]LHSImproved HS. It considers three main modifications: an adaptative global pitch adjustment, incorporation of an OBL technique, and a competition selection strategy.Ouyang et al. [142]SRHSSelective refining HS. The best harmony is considered as a representation of the current experience and knowledge of the HM. A tournament selection scheme is used. Instead of taking decisions based only on the value of the objective function, the value of every decision variable is considered.Shabani et al. [143]GOGHSGlobal HS with generalized OBL. Relevant information from the global-best harmony is extracted to create new candidate solutions. OBL is employed.Guo et al. [144]TPHSTwo-phase HS. The algorithm is divided in two phases, the first is highly focused in exploration by the implementation of the catastrophic mutation concept used in GA, and the second phase uses local search to improve exploitation.Assad and Deep [145]DSAHSEnhanced HS with dual strategies and adaptive parameters. An improved pitch adjustment operation with dual strategies is incorporated to enhance the exploitation ability of the algorithm. The operation may use a best harmony or a random harmony-guided search strategy according to a probability decision. Furthermore, the parameters of the algorithm are self-adjusted according to the new candidate fitness.Wang et al. [146]HS-SAInspired by the SA concept, a new parameter called temperature is used so that inferior harmonies could also be considered. Two additional SA parameters are employed.Assad and Deep [147]DHSDifferential-based HS. It considers the introduction of a position indicator to enable continuous HS to solve the discrete job shop scheduling problem. It also includes the incorporation of the DE mutation concept into the pitch adjustment step.Zhao et al. [148]CIHSChaotic improved HS. The concept of chaos is introduced into improved HS using logistic map to calculate PAR and bw. The best solutions are refined using a chaotic local search. New harmonies emulate the best harmony employing a PSO concept.Shafaati and Mojallali [149]ABHSAdaptative HS with best-based search strategy. It employs a best-based search strategy, which makes use of the global-best solution to enhance the search. HMCR and PAR are self-tuned based on information from the evolutionary search. The bw parameter is not considered in this algorithm.Guo et al. [150]DANGHSDynamic adjusting NGHS. Several dynamic parameter adjustment strategies to calculate the genetic mutation probability were used: linear, threshold, exponential, and concave-convex.Chiu et al. [151]AHSAmended HS. Two perturbation strategies are used: stochastic differential and OBL. The bw parameter is dynamically adjusted in terms of the mean value of a particular decision variable. A global dimension selection to balance exploration and exploitation is used.Ouyang et al. [152]MHS-PCLSParallel chaotic local search enhanced HS. It considers an initial point modification in PCLS to improve the performance of HS. It also uses MHS intersect mutation operation.Yi et al. [57]IHSGlobal learning and a modified random selection operation are proposed. A new improvement scheme based on the iteration modes (i.e., dimension-to-dimension, stochastic selection, vector, and matrix) is proposed. A total of four IHS variants were developed.Ouyang et al. [153]PHSβ-HC and imp. PHSβ-HCBest polynomial HS algorithm with bestβ-hill climbing. Two variants of the proposed method were introduced: PHSβ-HC and Imp. PHSβ-HC. The random selection step in HS is replaced by roulette wheel and tournament selection. A highly disruptive polynomial mutation replaces the pitch adjustment and random selection steps.Doush and Santos [154]EHS-CRPEnhanced HS with circular region perturbation. A global and local dimension selection is used by considering the global-best solution vector, and also local partial dimension variables are randomly selected. The selection operation is based on several PSO concepts. The circular region perturbation method is employed to improve search efficiency.Wu et al. [155]IMGHSAIntersect mutation global HS algorithm. Two methods were combined: intersect mutation and global HS. It considers the introduction of a new operator to reinforce exploration. HM is divided in two and a new operator considers both parts to generate new harmonies.Gholami et al. [156]SLHSSequence learning HS. It was adapted to solve the combinatorial flexible process planning problem. The optimal sequence is stablished by finding a proper successor for each operation.Luo [157]AHS-HCMAmeliorated HS with hybrid convergence mechanism. New harmonies are calculated through a convergence coefficient and concepts of PSO and DE based on the number of iterations. Non-linear dynamic domain in terms of current and maximum number of iterations when 1-HMCR.Zhu and Tang [158]MHSAModified HS algorithm. It considers the HM initialization by means of opposition numbers. Two new parameters inspired by the GWO method were introduced to calculate new harmonies. HMCR is calculated based on a Gaussian distribution. PAR is formulated as a non-linear exponential function.Gupta [159]AHS-DE-OBLAdaptative HS using DE and OBL. It includes the generation of new harmonies with OBL approach and the dynamic adjustment of HMCR, PAR, and search domain according to current and maximum number of iterations. The best and worst harmonies are used within a DE concept to calculate bw.Kang et al. [160]∗Variant of HS.Table 5
Harmony search variants based on hybridization with other metaheuristics.
AlgorithmRemarksReferenceSA-HSNelder–Mead simplex algorithm HS. The simplex method is incorporated in HM update using then best HM vectors + 1. The n best HM harmonies are kept in memory and the next vector is updated with the simplex solution. The HS improvisation stage considers these vectors along with the remaining harmonies in HM.Jang et al. [161]HHSAHybrid HS algorithm. Sequential quadratic programming (SQP) is used to enhance the local search performance. Two strategies are proposed: use of best harmonies in HM as initial points for the local search; performing a local search for each improvised solution vector. The strategy selection is based on a probabilistic parameter.Fesanghary et al. [162]HTHSAHybrid Taguchi HS algorithm. Two optimization phases are considered: (1) the Taguchi method is used to fine-tune decision variables intervals and to generate initial solutions and (2) HS is used to perform the global search.Yildiz [163]HS-DEDE is used to fine-tune HM elements considering all solution vectors as the individuals of the population.Gao et al. [164]HS-CSAThe clonal selection algorithm is used to improve the elements of the HM.Wang et al. [165]HS-solverIntegration of HS with non-linear excel solver algorithm. Two optimization approaches were proposed. The first one consists of the use of the local method once the HS solution is obtained. In the second approach, the algorithms operate simultaneously and new HS solutions are optimized with the local method at each iteration.Ayvaz et al. [166]HMGHSHybrid modified GHS. It was proposed to solve the block permutation flow shop scheduling problem. Continuous variables were coded to binary using the largest position value rule. HM initialization is based on problem’s heuristic rule. It includes a local search feature.Wang et al. [167]DMS-PSO-HSHS is incorporated into each subswarm of the dynamic multi-swarm particle swarm optimizer algorithm [168], in order to benefit from the particular exploration capabilities of the PSO technique and the exploitation abilities of HS.Zhao et al. [169]DCHSSADiscrete chaotic HS-based simulated annealing. New solutions of the SA algorithm are generated using HS parameters (HMCR and PAR). Chaotic mapping is used to generate new solutions when the probability is within HMCR.Askarzadeh [170]HS/BAHybridization of HS with the bat algorithm (BA) [171] where the mutation operator is modified by applying the concept of pitch adjustment to improve convergence.Wang and Guo [172]HS/FAHS and firefly algorithm (FA) [173]. Hybrid method that combines HS for exploration and FA for exploitation. It also incorporates HMCR and PAR to perform a pitch adjustment in FA mutation step.Guo et al. [174]AntHSHS and ACO. It considers the addition of two ACO concepts to HS: pheromone values and heuristics.Amini and Ghaderi [175]GHSACOGHS with ACO. A probabilistic factor is added to switch between using GHS or ACO for the improvisation step. The use of ACO was proposed to enhance exploitation.Fouad et al. [176]CPSO-MHSCooperative particle swarm optimizer-modified HS. It is a hybrid algorithm comprising a PSO variant (CPSO) and a modified HS version. MHS is employed for exploration and CPSO for exploitation. Note that PAR and bw parameters were removed using CPSO for local search.Zhang and Li [177]HS/CSHS with CS. It applies HS strategy in CS in order to increase diversity of cuckoo populations in the optimization process. Specifically, the HS improvisation strategy is used in the mutation stage of CS.Wang et al. [178]HSTLBOHS with teaching learning-based optimization [179]. A dynamic selection strategy is used to determine which algorithm will perform the search for a given iteration is used. It also includes a dynamic change of PAR and bw. In the improvisation step, a dynamic strategy considering elements of the worst harmony is implemented.Tuo et al. [180]MBHSMine blast HS. It employs MBA and HS for the exploration and exploitation phases, respectively. HS uses concepts from IHS [78] and PSF-HS to calculate PAR and HMCR in terms of iterations. A new iteration-dependent dynamic harmonic memory is used.Sadollah et al. [181]HS-GWOHS with grey wolf optimizer. GWO is used to control PAR and bw parameters during the search, and OBL enhances the improvisation step.Alomoush et al. [182]H-TSHS with TS. The neighborhood search capabilities of TS are exploited after HS improvisation step. It was applied to the dynamic parallel row ordering problem.Gong et al. [183]HHSDEHybrid harmony DE algorithm. It improved mutation operation in DE making use of HS improvisation step. HMS is dynamically updated according to current iteration. Use of a self-parameter adjustment strategy to control method’s parameters.Fu et al. [184]HABC-HSA-MAHybrid ABC and HS algorithm-based metaheuristic approach. These algorithms are used sequentially, and the final output of HS optimization process is the input for the ABC method.Jayalakshmi et al. [185]HS∗HS combined with L-BFGS. After generating HS initial population, a local search using the quasi-Newton method is performed for each harmony, keeping record of the best value in HM.Sorokovikov and Gornov [186]UKF-HSUnscented Kalman filter HS. It was developed to solve the finite-element model updating using single and multi-objective optimization.Naranjo-Pérez et al. [187]AHHSAdaptive-hybrid HS. The Java algorithm is used in HS improvisation step. HMCR and bw are randomly defined at the beginning of the algorithm.Yücel et al. [115]HSCAHSHybrid sine-cosine approach HS. HS is used in parallel with SCA to improve solution quality. HS generates a new population seeking to improve the worst harmonic in HM.Singh and Kaur [188]RDEGA-HSReward population-based differential genetic HS. It combines two algorithms proposed by the authors. The first, reward-population HS (RHS), works with a set of independent subharmonies and a reward harmony. The second, differential genetic HS (DEGA-HS), removes the bw parameter and employs two DE mutation strategies and a concept of the parthenogenetic algorithm in HS improvisation stage.Zhang et al. [189]∗Variant of HS.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 as(19)Optimizef1x,f2x,…,fmx,where 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.


## 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 via(8)PARiter=PARmin+PARmax−PARminitermax∗iter,(9)bwiter=bwmaxexplnbwmin/bwmaxitermax∗iter,where 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 (itermax), 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 xinew 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.
Step 2: HS algorithmStep 2: GHS algorithm(1)for each i∈1,ndo(2) ifrand1≤HMCRthen(3) xinew=xij where j∈1,HMS(4) ifrand2≤PARthen(5) xinew=xinew±rand∗bw(6) end if(7) else(8) xinew=randLBi,UBi(9) end if(10) end forfor each i∈1,ndoifrand1≤HMCRthenxinew=xij where j∈1,HMSifrand2≤PARthenxinew=xkbest where k∈1,nend ifelsexinew=randLBi,UBiend ifend forSuperscript best corresponds to the index of the best harmony in HM.


## 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 pm to avoid local solutions caused by premature convergence. Also, xnew replaces xHMS 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.
Step 2: NGHS algorithm(1)for each i∈1,ndo(2)  xR=2∗xibest−xiworst(3)  ifxR>UBithen(4)   xR=UBi(5)  else ifxR<LBithen(6)   xR=LBi(7)  end if(8)  xinew=xiworst+rand∗xR−xiworst(9)  ifrand≤pmthen(10)   xinew=randLBi,UBi(11)  end if(12) end forxR is the value of the decision variable with adaptative step. Superscript worst corresponds to the index of the worst harmony in HM.


## 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:(10)bwiter=bwmax−bwmax−bwminitermax∗2⋅iter,if iter<itermax2,bwmin,if iter≥itermax2.


## 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:(11)xinew=xij∈xi1,xi2,…,xiHMS+G0,1∗bw,with probability HMCR,randLBi,UBi,with probability1−HMCR,where a Gaussian distribution with mean 0 and standard deviation of 1 is used to calculate the random number G0,1. On the other hand, the pitch adjustment mechanism is given by(12)xinew=xibest+ϕ∗bw,with probability PAR∗HMCR,xinew,with probability1−PAR∗HMCR,where ϕ∈ [−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 xinew takes the value of xibest as in GHS. HMCR and PAR are dynamically adapted according to(13)HMCRiter=HMCRmin+HMCRmax−HMCRmin∗iteritermax,PARiter=PARmin+PARmax−PARmin∗iteritermax.


## 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 by(14)HMCRiter∈randn1+n,11+n,PARiter=1−iter−1itermax,bwi=xih−xim,ifxjh≠xji,UBi−LBi∗exp−UBi−LBi∗iteritermax,otherwise,where superscripts h and m refer to random integers within [1, HMS]. The pitch adjustment mechanism of ESGHS is given by(15)xinew=xih+ϕ∗xim−xih.Ifxim=xih at any point, the following equation is used for the pitch adjustment:(16)xinew=xih+ϕ∗UBi−LBi∗exp−UBi−LBi∗iteritermax.In the random selection (1 − HMCR), the new harmony element is obtained as follows:(17)xinew=randGμi,σ,where a Gaussian distribution with mean μi=xibest and standard deviation σ=1−iter−1/itermax2 is used to calculate the random number randG.


## 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 APiter is calculated at each iteration to decide if the new solution is accepted or rejected following the next criterion:(18)APiter=1,iffxnew≤fxHMS,fxHMS−fxbestfxnew−fxbest,iffxnew>fxHMS.


## 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 Tables3–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.
AlgorithmRemarksReferenceHS∗Based on a detailed mathematical analysis of the HM variance, the bw parameter was defined as the population standard deviation.Mukhopadhyay et al. [89]HS∗Self-adaptive HS. The bw parameter is replaced into the new harmony calculations considering the minimum and maximum values in HM. Furthermore, low-discrepancy sequences are used to initialize HM in place of random number generation.Wang and Huang [90]HS∗Adaptive HS. Dynamical adjustment of PAR and HMCR by a probabilistic selection based on mean values of the HM matrix.Hasançebi et al. [91]PSF-HSParameter-setting-free HS. Three steps were proposed to avoid the tuning of HS parameters: (1) a random tuning where HM is filled with random solution vectors, (2) a rehearsal step where a certain amount of solution vectors is generated and an additional HMS size memory keeps track of harmony source, and (3) a performance step where HMCR and PAR are used for each decision variable.Geem and Sim [92]PHSProposed HS. The bw parameter is calculated using derivatives (finite difference method) and lower and upper variable bounds. The probability of pitch adjustment is modified. An improvement to PHS was also presented where two bw definitions are used in the improvisation step. New decision probabilities are employed.Jaberipour and Khorram [93]EHSExplorative HS. The standard deviation of the current population defines the value of bw. Penalty strategy is used for constraint handling.Das et al. [94]NGHSA NGHS adaption to solve large-scale knapsack problems. The removal of HMCR and PAR is proposed to introduce a genetic mutation concept. A strategy from PSO is used to guide the worst solutions to global-best surroundings.Zou et al. [95]DHSDiscrete HS. It is designed to handle discrete variables. PAR and bw are set considering the defined bounds as well as the current and maximum number of iterations.Askarzadeh [96]NSHSNovel self-adaptive HS. HMCR and bw are fixed according to the problem dimension and decision variable bounds.Luo [97]LAHSLearning automata-based HS. Smart parameter setting of HMCR, PAR, and bw through a learning automata implementation [98].Enayatifar et al. [99]DSHSDynamic self-adaptive HS. It is based on the implementation of the best-to-worst ratio (BtW) concept in HM. PAR is dynamically adjusted in terms of BtW. Two bandwidth parameters are used, active bw and dynamic bandwidth value dbw. These parameters are calculated for each decision variable and are dynamically adjusted according to its standard deviation in all the HM.Kattan and Abdullah [100]ABHSAdaptative binary HS. New harmonies are created based on single or multiple vectors in HM. PAR is modified according to the best HM vector. HMCR and PAR are obtained using several mechanisms: linear decrease, linear and non-linear increase, and random, where all are calculated as a function of current and total iterations. ABHS can generate one or more harmonies at each generation. HM update is performed in serial or parallel according to the number of generated new candidates.Wang et al. [101]ABHSAdjustable bandwidth HS. Minimum and maximum values for PAR and HMCR are defined. Additionally, a saturation parameter (SP) is proposed to control bw parameter. An alternative version of ABHS with changes in bw calculation is also presented.Contreras et al. [102]GDHSGlobal dynamic HS. It is based on a dynamic change in bw, PAR, and HMCR parameters, as well as in the domain of HM variables.Khalili et al. [103]PAHSParameter adaptive HS. Four dynamic change strategies for PAR and HMCR in terms of the number of iterations were implemented.Kumar et al. [104]SFHSSelf-regulated fretwidth HS. It considers changes in bw (also called fretwidth parameter) at each iteration based on several criteria.Amaya et al. [105]SBHSSimplified binary HS. It focuses on changes to the improvisation step. The method is parameter free (i.e., PAR and bw were removed), and HMCR is dynamically adjusted. It incorporates greedy local search. SBHS was used to solve large-scale knapsack problems.Kong et al. [106]HS∗Dynamic adjustment of bw. For a particular decision variable, bw is defined according to the current iteration and its maximum and minimum values.Kalivarapu et al. [107]GGHSGaussian GHS. HM is initialized based on a Gaussian distribution and it is adjusted considering a dynamic bw. The best harmony is also adjusted using a dynamic standard deviation.Keshtegar and Sadeq [108]ImHSImproved HS. The bw parameter decrements according to an exponential function in terms of the number of iterations. A selection scheme based on a proposed parameter called rate of improvement adjustment is applied.Portilla-Flores et al. [109]SSaHSSemi-self-adaptive HS. The bw parameter is updated according to a self-adaptative rule, which considers the information about the current population. Two new parameters are stablished by an empirical analysis.Zhao et al. [110]PSFHSModified parameter-setting-free HS. The number of harmonies in HM is modified by a positive integer coefficient. PAR and HMCR domains are redefined considering the minimum and maximum values, respectively. Recommended limits for HMCR and PAR were provided.Shaqfa and Orbán [111]BWM_HSBest-worst-mean HS. It considers the replacement of HMCR for three novel criteria to enhance exploration/exploitation balance. A self-adaptative probabilistic factor is incorporated to choose between rules.Talaei et al. [112]FHSFuzzy HS. A fuzzy system is used to handle HMCR and PAR throughout iterations.Peraza et al. [113]EHSEnhanced HS. It considers dynamic adjustment of PAR and bw. PSO strategy is implemented in new harmony generation.Niu et al. [114]AHSAdaptative HS. Parameters HMCR and bw are defined as function of current and total number of iterations.Yücel et al. [115]SIHSSelf-adaptative improved HS. It considers the dynamic modification of bw starting with a large value and adjusting it at each generation.Houari et al. [116]∗Variant of HS.Table 4
Harmony search variants based on the algorithm structure modification.
AlgorithmRemarksReferenceHS∗Algorithm for structural optimization. It implies the initialization of HM memory using the finite element method considering problem’s constraints. This approach is also used during the improvisation step to assess new harmonies.Lee and Geem [117]HS∗It was proposed to solve discrete variable optimization problems employing a stochastic derivative formulation based on several probability criteria using HS parameters. It was tested on benchmark and real-world problems.Geem [118]HS∗HS with SA concept. It is based on the modification of PAR and bw. Starting from a maximum initial value, the parameters are slowly reduced imitating the cooling concept of the SA algorithm.Taherinejad [119]IHSImproved HS. It considers the dynamic modification of PAR based on dispersed particle swarm optimization (DPSO) [120]. The information index, which relates the difference of the particles according to its fitness value, is incorporated.Coelho and Bernert [121]DHSDifferential HS. DE mutation strategy is implemented in the improvisation step.Chakraborty et al. [122]PSHSParticle swarm HS. A concept of particle swarm optimization is used in the improvisation step.Geem [123]CHSIt uses chaotic maps (e.g., logistic, tent, Gauss, and so on) to substitute generally employed random numbers in HS operators. The algorithm was tested in multi-modal benchmark functions.Alatas [124]GHS + LEMMachine learning tools are employed to create new populations in HS according to the learnable evolution model (LEM) [125]. Also, an evolutionary related method called the Darwinian is used to decide which members of the population are better to fulfil certain tasks.Cobos et al. [126]ITHSIntelligent tuned HS. HM is divided in two groups in order to balance exploration and exploitation. The first one contains those solutions with an objective value less or equal to the mean value in HM, and the second group includes the rest of solution vectors. PAR is dynamically adjusted.Yadav et al. [127]NDHSA tournament selection rule is used to choose the candidate harmony from the HM. The solution vectors with a better objective function value will be more likely to be used to generate new solutions. Also, the PAR and bw parameters are dynamically adjusted.Chen et al. [128]HSTLHS based on teaching-learning. It was proposed to solve complex high dimension optimization problems. It considers the introduction of a teaching-learning concept into the improvisation step and a dynamical adjustment of all parameters (i.e., HMCR, PAR, bw, and TLP) in terms of generations.Tuo et al. [129]IGHSImproved GHS. The opposition-based learning (OBL) [130] was used as an alternative initialization procedure to improve the quality of the initial memory vectors. In addition, a mutation strategy was implemented in place of the GHS pitch adjustment. The random selection step was substituted for a modification used in artificial bee colony (ABC) algorithm [131].Xiang et al. [132]IGHSIntelligent global HS. In the improvisation step, the new harmony mimics an element of the best solution vector stored in HM.Valian et al. [133]GSHSGeometric selective HS. Two harmonies are combined to create theith element of the new harmony using a tournament selection strategy. A mutation operator is implemented using PAR.Castelli et al. [134]OHSOpposition-based HS. OBL strategy is incorporated in the HM initialization and update steps.Banerjee et al. [135]HS∗Several mutation strategies are incorporated in the random consideration process (i.e., random, boundary, non-uniform, MPT, power, and polynomial).Hasan et al. [136]QOHSQuasi-oppositional HS. It considers the introduction of the quasi-oppositional-based learning strategy in HS and the application of the quasi-opposition concept for HM initialization. PAR and bw are obtained according to current generation. It also implies the improvisation of new harmonies and populations on quasi-oppositional basis and jumping probability according to search space bounds.Shiva et al. [137]iHSIsland HS. It considers the incorporation of the island model [138] in the HM memory.Al-Betar et al. [139]DH/bestRandom initialization of HM is replaced by a new technique, which reduces randomness. The conventional pitch adjustment method is substituted by a DE mutation strategy.Abedinpourshotorban et al. [140]MHSModified HS. The bw parameter is removed. An intersect operation considering PAR is included, and after each improvisation, the algorithm performs a cellular local search.Yi et al. [141]LHSImproved HS. It considers three main modifications: an adaptative global pitch adjustment, incorporation of an OBL technique, and a competition selection strategy.Ouyang et al. [142]SRHSSelective refining HS. The best harmony is considered as a representation of the current experience and knowledge of the HM. A tournament selection scheme is used. Instead of taking decisions based only on the value of the objective function, the value of every decision variable is considered.Shabani et al. [143]GOGHSGlobal HS with generalized OBL. Relevant information from the global-best harmony is extracted to create new candidate solutions. OBL is employed.Guo et al. [144]TPHSTwo-phase HS. The algorithm is divided in two phases, the first is highly focused in exploration by the implementation of the catastrophic mutation concept used in GA, and the second phase uses local search to improve exploitation.Assad and Deep [145]DSAHSEnhanced HS with dual strategies and adaptive parameters. An improved pitch adjustment operation with dual strategies is incorporated to enhance the exploitation ability of the algorithm. The operation may use a best harmony or a random harmony-guided search strategy according to a probability decision. Furthermore, the parameters of the algorithm are self-adjusted according to the new candidate fitness.Wang et al. [146]HS-SAInspired by the SA concept, a new parameter called temperature is used so that inferior harmonies could also be considered. Two additional SA parameters are employed.Assad and Deep [147]DHSDifferential-based HS. It considers the introduction of a position indicator to enable continuous HS to solve the discrete job shop scheduling problem. It also includes the incorporation of the DE mutation concept into the pitch adjustment step.Zhao et al. [148]CIHSChaotic improved HS. The concept of chaos is introduced into improved HS using logistic map to calculate PAR and bw. The best solutions are refined using a chaotic local search. New harmonies emulate the best harmony employing a PSO concept.Shafaati and Mojallali [149]ABHSAdaptative HS with best-based search strategy. It employs a best-based search strategy, which makes use of the global-best solution to enhance the search. HMCR and PAR are self-tuned based on information from the evolutionary search. The bw parameter is not considered in this algorithm.Guo et al. [150]DANGHSDynamic adjusting NGHS. Several dynamic parameter adjustment strategies to calculate the genetic mutation probability were used: linear, threshold, exponential, and concave-convex.Chiu et al. [151]AHSAmended HS. Two perturbation strategies are used: stochastic differential and OBL. The bw parameter is dynamically adjusted in terms of the mean value of a particular decision variable. A global dimension selection to balance exploration and exploitation is used.Ouyang et al. [152]MHS-PCLSParallel chaotic local search enhanced HS. It considers an initial point modification in PCLS to improve the performance of HS. It also uses MHS intersect mutation operation.Yi et al. [57]IHSGlobal learning and a modified random selection operation are proposed. A new improvement scheme based on the iteration modes (i.e., dimension-to-dimension, stochastic selection, vector, and matrix) is proposed. A total of four IHS variants were developed.Ouyang et al. [153]PHSβ-HC and imp. PHSβ-HCBest polynomial HS algorithm with bestβ-hill climbing. Two variants of the proposed method were introduced: PHSβ-HC and Imp. PHSβ-HC. The random selection step in HS is replaced by roulette wheel and tournament selection. A highly disruptive polynomial mutation replaces the pitch adjustment and random selection steps.Doush and Santos [154]EHS-CRPEnhanced HS with circular region perturbation. A global and local dimension selection is used by considering the global-best solution vector, and also local partial dimension variables are randomly selected. The selection operation is based on several PSO concepts. The circular region perturbation method is employed to improve search efficiency.Wu et al. [155]IMGHSAIntersect mutation global HS algorithm. Two methods were combined: intersect mutation and global HS. It considers the introduction of a new operator to reinforce exploration. HM is divided in two and a new operator considers both parts to generate new harmonies.Gholami et al. [156]SLHSSequence learning HS. It was adapted to solve the combinatorial flexible process planning problem. The optimal sequence is stablished by finding a proper successor for each operation.Luo [157]AHS-HCMAmeliorated HS with hybrid convergence mechanism. New harmonies are calculated through a convergence coefficient and concepts of PSO and DE based on the number of iterations. Non-linear dynamic domain in terms of current and maximum number of iterations when 1-HMCR.Zhu and Tang [158]MHSAModified HS algorithm. It considers the HM initialization by means of opposition numbers. Two new parameters inspired by the GWO method were introduced to calculate new harmonies. HMCR is calculated based on a Gaussian distribution. PAR is formulated as a non-linear exponential function.Gupta [159]AHS-DE-OBLAdaptative HS using DE and OBL. It includes the generation of new harmonies with OBL approach and the dynamic adjustment of HMCR, PAR, and search domain according to current and maximum number of iterations. The best and worst harmonies are used within a DE concept to calculate bw.Kang et al. [160]∗Variant of HS.Table 5
Harmony search variants based on hybridization with other metaheuristics.
AlgorithmRemarksReferenceSA-HSNelder–Mead simplex algorithm HS. The simplex method is incorporated in HM update using then best HM vectors + 1. The n best HM harmonies are kept in memory and the next vector is updated with the simplex solution. The HS improvisation stage considers these vectors along with the remaining harmonies in HM.Jang et al. [161]HHSAHybrid HS algorithm. Sequential quadratic programming (SQP) is used to enhance the local search performance. Two strategies are proposed: use of best harmonies in HM as initial points for the local search; performing a local search for each improvised solution vector. The strategy selection is based on a probabilistic parameter.Fesanghary et al. [162]HTHSAHybrid Taguchi HS algorithm. Two optimization phases are considered: (1) the Taguchi method is used to fine-tune decision variables intervals and to generate initial solutions and (2) HS is used to perform the global search.Yildiz [163]HS-DEDE is used to fine-tune HM elements considering all solution vectors as the individuals of the population.Gao et al. [164]HS-CSAThe clonal selection algorithm is used to improve the elements of the HM.Wang et al. [165]HS-solverIntegration of HS with non-linear excel solver algorithm. Two optimization approaches were proposed. The first one consists of the use of the local method once the HS solution is obtained. In the second approach, the algorithms operate simultaneously and new HS solutions are optimized with the local method at each iteration.Ayvaz et al. [166]HMGHSHybrid modified GHS. It was proposed to solve the block permutation flow shop scheduling problem. Continuous variables were coded to binary using the largest position value rule. HM initialization is based on problem’s heuristic rule. It includes a local search feature.Wang et al. [167]DMS-PSO-HSHS is incorporated into each subswarm of the dynamic multi-swarm particle swarm optimizer algorithm [168], in order to benefit from the particular exploration capabilities of the PSO technique and the exploitation abilities of HS.Zhao et al. [169]DCHSSADiscrete chaotic HS-based simulated annealing. New solutions of the SA algorithm are generated using HS parameters (HMCR and PAR). Chaotic mapping is used to generate new solutions when the probability is within HMCR.Askarzadeh [170]HS/BAHybridization of HS with the bat algorithm (BA) [171] where the mutation operator is modified by applying the concept of pitch adjustment to improve convergence.Wang and Guo [172]HS/FAHS and firefly algorithm (FA) [173]. Hybrid method that combines HS for exploration and FA for exploitation. It also incorporates HMCR and PAR to perform a pitch adjustment in FA mutation step.Guo et al. [174]AntHSHS and ACO. It considers the addition of two ACO concepts to HS: pheromone values and heuristics.Amini and Ghaderi [175]GHSACOGHS with ACO. A probabilistic factor is added to switch between using GHS or ACO for the improvisation step. The use of ACO was proposed to enhance exploitation.Fouad et al. [176]CPSO-MHSCooperative particle swarm optimizer-modified HS. It is a hybrid algorithm comprising a PSO variant (CPSO) and a modified HS version. MHS is employed for exploration and CPSO for exploitation. Note that PAR and bw parameters were removed using CPSO for local search.Zhang and Li [177]HS/CSHS with CS. It applies HS strategy in CS in order to increase diversity of cuckoo populations in the optimization process. Specifically, the HS improvisation strategy is used in the mutation stage of CS.Wang et al. [178]HSTLBOHS with teaching learning-based optimization [179]. A dynamic selection strategy is used to determine which algorithm will perform the search for a given iteration is used. It also includes a dynamic change of PAR and bw. In the improvisation step, a dynamic strategy considering elements of the worst harmony is implemented.Tuo et al. [180]MBHSMine blast HS. It employs MBA and HS for the exploration and exploitation phases, respectively. HS uses concepts from IHS [78] and PSF-HS to calculate PAR and HMCR in terms of iterations. A new iteration-dependent dynamic harmonic memory is used.Sadollah et al. [181]HS-GWOHS with grey wolf optimizer. GWO is used to control PAR and bw parameters during the search, and OBL enhances the improvisation step.Alomoush et al. [182]H-TSHS with TS. The neighborhood search capabilities of TS are exploited after HS improvisation step. It was applied to the dynamic parallel row ordering problem.Gong et al. [183]HHSDEHybrid harmony DE algorithm. It improved mutation operation in DE making use of HS improvisation step. HMS is dynamically updated according to current iteration. Use of a self-parameter adjustment strategy to control method’s parameters.Fu et al. [184]HABC-HSA-MAHybrid ABC and HS algorithm-based metaheuristic approach. These algorithms are used sequentially, and the final output of HS optimization process is the input for the ABC method.Jayalakshmi et al. [185]HS∗HS combined with L-BFGS. After generating HS initial population, a local search using the quasi-Newton method is performed for each harmony, keeping record of the best value in HM.Sorokovikov and Gornov [186]UKF-HSUnscented Kalman filter HS. It was developed to solve the finite-element model updating using single and multi-objective optimization.Naranjo-Pérez et al. [187]AHHSAdaptive-hybrid HS. The Java algorithm is used in HS improvisation step. HMCR and bw are randomly defined at the beginning of the algorithm.Yücel et al. [115]HSCAHSHybrid sine-cosine approach HS. HS is used in parallel with SCA to improve solution quality. HS generates a new population seeking to improve the worst harmonic in HM.Singh and Kaur [188]RDEGA-HSReward population-based differential genetic HS. It combines two algorithms proposed by the authors. The first, reward-population HS (RHS), works with a set of independent subharmonies and a reward harmony. The second, differential genetic HS (DEGA-HS), removes the bw parameter and employs two DE mutation strategies and a concept of the parthenogenetic algorithm in HS improvisation stage.Zhang et al. [189]∗Variant of HS.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 as(19)Optimizef1x,f2x,…,fmx,where 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, Figure5 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.
AlgorithmApplicationReferenceHSCombined heat and power economic dispatch problemVasebi et al. [214]HSDesign of shell and tube heat exchangersFesanghary [73]IHSAir cooled heat exchanger designDoodman et al. [215]HSWater distribution system designGeem [216]HS-SQPHeat exchanger network designKhorasany and Fesanghary [217]HSParameter estimation for phase equilibrium modelingBonilla-Petriciolet et al. [218]NGHSChemical equation balance problemZou et al. [219]HS, IHS, GHSPhase stability analysis and phase equilibrium calculations in non-reactive systemsBonilla-Petriciolet et al. [220]HSDistillation sequence designCabrera-Ruíz et al. [221]HSParameter estimation for phase equilibrium modelingMerzougui et al. [222]HS, IHS, GHSPhase equilibrium calculations in reactive systemsBonilla-Petriciolet et al. [223]HS, IHS, GHSParameter estimation for phase equilibrium modelingBonilla-Petriciolet [224]HSLife cycle cost analysisFesanghary et al. [197]HSParameter estimation for the modeling of electrolyte solutionsBonilla-Petriciolet et al. [225]IHSDesign of plate-fin heat exchangersYousefi et al. [226]APF-HSGround water pollution problemJiang et al. [227]HSAEAReaction kinetic parameter estimationMa et al. [228]I-ITHSDesign of shell and tube heat exchangersTurgut et al. [229]HSLife cycle cost analysisAsadi [230]HSParameter estimation for the equilibrium modeling of biofuel-based systemsMerzougui et al. [231]HSWater distribution system designYoo et al. [232]HSParameter estimation for the phase equilibrium modelingMerzougui et al. [233]HSPhase stability analysis, reactive and non-reactive phase equilibrium calculations, and parameter estimationFernández-Vargas et al. [234]IHSHeat integrated distillation system design in olefin productionLashkajani et al. [235]DHSLife cycle cost analysisMaleki et al. [236]HSMOPump scheduling problemDe Paola et al. [237]HyHSWater distribution system designJung et al. [238]HSThermal cracking processBoto et al. [239]HSWater distribution system designYoo et al. [240]HSHeat integrated distillation system designSukpancharoen et al. [241]CHSDiesel blending problemGao et al. [242]SFHSParameter estimation for the phase equilibrium modeling of supercritical fluidsAmaya et al. [243]HSSeparation cascade distillation designMansourzadeh et al. [244]HSReactive azeotrope calculation and prediction of dew point pressures in binary double retrograde vaporization systemsPlatt et al. [245]HS-GWOPhase stability analysisHernández-Pérez et al. [246]HSParameter estimation for phase equilibrium modelingReggab et al. [247]HSPhase stability analysisSmejkal et al. [248]HSPhase stability analysisFonseca-Pérez et al. [249]HSMulti-effect evaporator designYadav et al. [250]HS-GSInverse chemical kinetic problemEnikeeva et al. [251]Table 7
Problem description of thermodynamic calculations.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferenceParameter estimation for VLE modelingGlobal optimizationMinimize∑j=1ndat∑i=1cγijexp−γijcalc/γijexp2, where γij is the activity coefficient, c is the number of components, ndat is the number of VLE data, exp is the experimental value, and calc is the calculated value.Binary interaction parameters (Aij) for local composition models. (a) Wilson model: −10,000≤Aij≤20,000. (b) NRTL model: −2,000≤Aij≤5,000. (c) UNIQUAC model: −5,000≤Aij≤20,000.Binary interaction parameters (Aij and αij for local composition models. (a) Wilson model: −8,500≤Aij≤320,000. (b) UNIQUAC model: −5,000≤Aij≤20,000. (c) NRTL model: −2,000≤Aij≤5,0000.01≤αij≤10.UnconstrainedUnconstrainedBonilla-Petriciolet et al. [218]Bonilla-Petriciolet [224]; Fernández-Vargas et al. [234]Phase stability analysisGlobal optimizationMinimize∑i=1cyiμi|y−μi|z, where yi is the mole fraction of trial phase, μi is the chemical potential, z is the feed composition, and c is the number of components.Minimize ∑i=1nCiμiC1,…,Cc−μix∗−PC,…,Cc−Px∗, where Ci is the component concentration, μi is the chemical potential, c is the number of components, P is the pressure, and x∗ is the initial concentration.The constrained problem with decision variablesyi (mole fractions) was formulated as an unconstrained problem with the use of new decision variables 0≤βi≤1 related to compositions via the equations: niy=βizinF and yi=niy/∑j=1cnjy, where niy is the mole number of trial composition and nF is the total number of moles in the feed composition under analysis. The optimization problem dimension is given by c.Component concentrations Ci, Ci≥0,i=1,…,c, ∑i=1cbiCi<1, where bi is a PR equation parameter. The solution y=C,…,c was transformed with the next equation xi=1−εyi/bi into the feasible domain xi≥0, ∑i=1cxi≤1, where ε is a small positive number and x is the transformed solution.UnconstrainedUnconstrainedBonilla-Petriciolet et al. [220]; Fernández-Vargas et al. [234]; Hernández-Pérez et al. [246]; Fonseca-Pérez et al. [249]Smejkal et al. [248]Parameter estimation for LLE modelingGlobal optimizationMinimize∑k=1mtie∑j=12∑i=1cxikjcal−xikjexp2, where xijk is the liquid mole fraction, c is the number of components, mtie is the number of tie lines, cal is the calculated value, and exp is the experimental valueSix NRTL binary interaction parametersAij: −1,000≤Aij≤2,000NRTL and UNIQUAC binary interaction parameters Aij of biofuel systems: −1,000≤Aij≤2,000NRTL and UNIQUAC binary interaction parameters Aij for food industry systems. For ternary and quaternary systems, 5 and 9 decision variables were determined, respectively. The remaining interaction parameters were calculated with the following relationship: ∑i=1c∑j=1j≠icAij−Aji=0, −1,000≤Aij≤2,000NRTL binary interaction parameters Aij for a ternary mixtureUnconstrainedUnconstrainedUnconstrainedUnconstrainedMerzougui et al. [222]Merzougui et al. [231]Merzougui et al. [233]Reggab et al. [247]Phase equilibrium calculations in reactive systemsGlobal optimizationMinimizeΔg−∑j=1πlnKeqN−1nref,j, where Δg=∑j=1π∑i=1cnijlnxijγij, Δg is the Gibbs free energy of mixing, lnKeq is a row vector of logarithms of chemical equilibrium constants for all independent chemical reactions, N is an invertible square matrix containing the stoichiometric coefficients vi of a set of reference components, nref is a column vector containing the number of moles of each reference component, nij is the number of moles, xij is the mole fraction, γij is the activity coefficient, c is the number of components, and π is the number of phases.Number of molesnij. To reduce the number of decision variables, the following expression was used: niπ=niF−viN−1nref,F−nref,π − ∑j=1π−1nij−viN−1nref,j, i=1,…,c−r. The number of decision variables is given by cπ−1+r, where r is the number of independent chemical reactions. nij>0,i=1,…,c;j=1,…,π, where vi is the vector of the stoichiometric coefficients of component i.Mass balance restrictions:∑j=1πnij−viN−1nref,j=niF−viN−1nref,Fi=1,…,c−r, where niF is the number of moles in the feed.Bonilla-Petriciolet et al. [223]; Fernández-Vargas et al. [234]Parameter estimation for VLE modeling with error-in-variable approachGlobal optimizationMinimize∑j=1ndat∑i=1cxijt−xij2/σxi2+yijt−yij2/σyi2 + ∑j=1neTjt−Tj2/σT2+Pjt−Pj2/σP2, where xijt is the true value of the liquid composition, yijt is the true value of the vapor composition, Tjt is the true value of the temperature, Pjt is the true value of the pressure, xij is the experimental value of the liquid composition, yij is the experimental value of the vapor composition, Tj is the experimental value of the temperature, Pj is the experimental value of the pressure, c is the number of components, ndat is the number of experimental data, σxi is the standard deviation associated to the liquid composition, σyi is the standard deviation associated to the vapor composition, σT is the standard deviation associated to the temperature, and σP is the standard deviation associated to the pressure.The set of true values for each state variable (i.e.,xijt, yijt, Tjt, Pjt) and the Wilson model parameters Aij. The number of decision variables is given by npar+c∗ndat where npar is the number of model parameters. −10,000≤Aij≤200,000, xij−3σ≤xijt≤xij+3σ, T−3σ≤Tt≤T+3σ. The dimensionality of the VLE problem was reduced by the elimination of Pt and yit in the objective function according to P=∑i=1cγixiPi0, yi=γixiPi0/∑j=1cγjxjPj0,i=1,…,c where γi is the activity coefficient and Pi0 is the vapor pressure at system temperature.UnconstrainedBonilla-Petriciolet [224]Modeling activity coefficients for electrolyte solutionsGlobal optimizationMinimize∑i=1ndatlnγ±expi−lnγ±calci2, where γ± is the electrolyte mean activity coefficient, ndat is the number of experimental data, exp is the experimental value, and calc is the calculated value.eNRTL interaction parametersτcas and τsca: −50≤τcas≤10, −10≤τsca≤50.UnconstrainedBonilla-Petriciolet et al. [225]Gibbs free energy minimization in non-reactive systemsGlobal optimizationMinimize∑j=1π∑i=1cnijμij, where nij is the number of moles, μij is the chemical potential, c is the number of components, and π is the number of phases.Introduction of new decision variables0≤βi≤1 instead of nij to transform the problem from constrained to unconstrained according to the following expressions: ni1=βi1zinF, nij=βijzinF−∑m=1j−1nim, niπ=zinF−∑m=1π−1nim, i=1,…,c;j=2,…,π−1, where zi is the composition of the feed and nF is the number of moles in the feed. The number of decision variables is cπ−1.UnconstrainedFernández-Vargas et al. [234]Parameter estimation for the VLE modeling of supercritical fluidsGlobal optimizationMinimize∑indatPiexp−Picalc2, where ndat is the number of experimental data, Pi is the pressure, exp is the experimental value, and calc is the calculated value.PR and Van Laar parameters.UnconstrainedAmaya et al. [243]Reactive azeotrope calculationGlobal optimizationMinimize∑k=15Fk2, where Fk is the residue of the expressions, Pyi−xiγiPisat=0,i=1,2, P1−y1−y2−1−x1−x2γ3P3sat=0, logK−vtloga^1loga^2⋮loga^c=0, X1−Y1=0, P is the system pressure, yi is the vapor molar fraction, xi is the liquid molar fraction, γi is the activity coefficient, Pisat is the pure component saturation pressure, vt is a vector with the stoichiometric coefficients, K is the chemical equilibrium constant, a^i is the activity of component i, X1 represents the liquid transformed coordinates, and Y1 represents the vapor transformed coordinates.Vapor and liquid molar fractions,y1, y2, x1, x2, and system temperature T.K=∏i=1ca^ivi, ∑i=1cxi=1, ∑i=1cyi=1, where c is the number of components.Platt et al. [245]Dew point pressure prediction for double retrograde vaporization systemsGlobal optimizationMinimize∑i=12Fi2, F1=ϕ^1Lx1−ϕ^1Vy1=0, F2=ϕ^2L1−x1−ϕ^2Vy2=0, where Fi is a function of the non-linear system of equations, ϕ^i is the fugacity coefficient, xi is the liquid molar fraction, yi is the vapor molar fraction, L is the liquid phase, and V is the vapor phase.System pressureP and liquid molar fraction x1.∑i=1cxi=1, ∑i=1cyi=1, where c is the number of components.Platt et al. [245]Table 8
Problem description of the design of heat exchangers.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferenceShell and tube heat exchanger designGlobal optimizationMinimizeCTotal, CTotal=Coperating+Ccapital, Ccapital=CciR1+iRTL/1+iRTL−1, Coperating=TP⋅ec⋅EΔP, where CTotal is the annual total cost, Ccapital is the annual capital cost, Coperating is the annual operating cost, TP is the yearly period of operation time, EΔP is the pumping power, ec is the energy cost per unit, Cc is the capital cost, iR is the interest rate, and TL is the technical life.Discrete and continuous variables: inside shell diameter, tube outside diameter, tube pitch, length, baffle cut, central baffle spacing, number of tube passes, tube layout characteristic angle, and number of tubes.No information available about constraints.Fesanghary et al. [252]MinimizeCtot, Ctot=Ci+C0D, Ci=a1+a2Sa3, C0=PCeH, C0D=∑k=1nyC0/1+ik, where Ctot is the total cost of the heat exchanger, Ci is the capital investment, C0D is the discounted operation cost, C0 is the operation cost, Ce is the energy cost, a1−3 represent the heat exchanger coefficients, S is the heat exchanger surface area, P is the efficiency factor, H is the annual number of work hours, ny is the number of years, and i is the annual discount rate.Discrete and continuous variables: the number of tube passes, shell and tube outside diameter, and baffle spacing.ΔPt=ΔPmajor+ΔPminor where ΔPt is the tube side pressure drop, ΔPmajor is the major loss pressure drop, and ΔPminor is the minor loss pressure drop.Turgut et al. [229]Heat exchanger network designGlobal optimizationMinimize∑i∈HPCCUQcu,i+∑j∈CPCHUQhu,j + ∑k∈NEC0+C1Akm+∑j∈CPC2+C3Ajm + ∑i∈HPC2+C3Aim. The first two terms correspond to the annual utility cost, while the rest of the terms represent the capital cost. CCU is the cold utility cost. CHU is the hot utility cost. Qcu,i is the heat exchanged between stream i and cold utility. Qhu,j is the heat exchanged between stream j and hot utility. C1−3 represent equipment dependent constants. A is the area. HP is the index set of hot process stream. CP is the index set of cold process stream. NE is the index set of process heat exchangers.The stream split fractions and the heat load of the exchangers.(1) Heat balance:Tiin−TioutFCpi = ∑k∈NL∑j∈CPQijk+Qcu,i, Tiout−TiinFCpj = ∑k∈NL∑i∈HPQijk+Qhu,j, where Ti is the temperature, in is the inlet value out is the outlet value, FCpi is the heat capacity flow rate, Qijk is the heat exchanged, NL is the index set of exchanger levels. (2) Non-negative bounds: Qijk≥0, Qcu,i≥0, Qhu,j≥0, y≥0, where y is the mole fraction vector. (3) Minimum temperature: dtij≥ΔTmin, dtcui≥ΔTmin, dthuj≥ΔTmin, where dt is the temperature approach and ΔTmin is the minimum temperature approach. (4) Mass conservation: ∑z∈SSyz=1, where SS is the index set of stream splits.Khorasany and Fesanghary [217]Air cooled heat exchanger designGlobal optimizationMinimizeCTotal, CTotal=Coperating+Ccapital, Coperating=TP⋅ec⋅Pfan, Ccapital=CciR1+iRTL/1+iRTL−1, Cc=aAb, where CTotal is the annual total cost, Ccapital is the annual capital cost, Coperating is the annual operating cost, TP is the yearly period of operation time, ec is the energy cost per unit, Pfan is the fan electrical power consumption, Cc is the air cooled heat exchanger capital cost, a,b are the cost factors, A is the tube area, iR is the interest rate, and TL is the technical life.Discrete and continuous variables: air velocity, tube material, length and outside diameter, fin height, number of tubes per row, and tube passes.No information available about constraints.Doodman et al. [215]Plate-fin heat exchanger designGlobal optimizationMinimizeAHT, AHT=Aa+Ab, Aa=LaLbNa1+2naHa−ta, Ab=LaLbNa1+2nbHb−tb, where AHT is the heat transfer area, La is the length of side a, Lb is the length of side b, Na is the number of layers in the hot side, H is the fin height, t is the fin thickness, and n is the fin frequency.Discrete and continuous variables:La, Lb, Na, H, t, n, and the lance length lfa.Constraints were handled using a penalty method.Q=εCminTa,1−Tb,1, ΔPa=2faLaGa2/ρaDh,a, ΔPb=2fbLbGb2/ρbDh,b, Lc=1.5, where Q is the heat transfer rate, ε is the effectiveness, Cmin is the minimum heat capacity rate, T is the temperature, f is the friction factor, G is the mass flow velocity, ρ is the density, Dh is the hydraulic diameter, Lc is the non-flow length, and a,b are fluids.Yousefi et al. [226]MinimizeΔPa/Pa,max+ΔPb/Pb,max, where ΔP is the pressure drop by friction, P is the pressure, a,b are fluids, and max is the maximum value.Discrete and continuous variables: the length of sidesa and b, number of layers in the hot side, fin height, thickness and frequency, and lance length.Constraints were handled using a penalty method.Q=εCminTa,1−Tb,1, ΔPa=2faLaGa2/ρaDh,a, ΔPb=2fbLbGb2/ρbDh,b, Lc=1.5, where Q is the heat transfer rate, ε is the effectiveness, Cmin is the minimum heat capacity rate, T is the temperature, f is the friction factor, L is the length, G is the mass flow velocity, ρ is the density, Dh is the hydraulic diameter, and Lc is the non-flow length.Yousefi et al. [226]Table 9
Problem description of distillation design and optimization.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferenceDesign of three distillation sequencesGlobal optimizationFor multi-component distillation: minimizeQ=fR,P0,TF,NF,NT, where Q is the column heat duty, R is the reflux ratio, P0 is the pressure of the column, TF is the temperature of the feed, NF is the feed stage number, and NT is the number of column stages.Discrete and continuous variables:R,P0,TF,NF, and NT.ym≥xm, where ym is the resulting purity vector and xm is the required purity vector.Cabrera-Ruiz et al. [221]For thermally coupled reactive distillation: minimizeQ=fR,P0,FD,FS,FL,NF1,NF2,N0Re,NfRe,NFL,NFV,NT1,NT2, where R is the reflux ratio, P0 is the pressure of the column, FD is the distillate flux, FS is the side flux, FL is the interconnection flow, NF1−2 represent the numbers of the feed stages, N0Re is the first reaction tray location, NfRe is the last reaction tray location, NFL is the location of the interconnection flow, NFV is the location of the vapor stream, NT1 is the number of stages of the main column, and NT2 is the number of stages of the stripper.Discrete and continuous variables:R,P0,FD,FS,FL,NF1,NF2,N0Re,NfRe,NFL,NFV,NT1, and NT2.ym≥xm, where ym is the resulting purity vector and xm is the required purity vector.Cabrera-Ruiz et al. [221]For dividing wall distillation: minimizeQ=fR,P0,FD,FS1,FS2,FL1,FV2,NF,NP0,NP,NS1,NS2,NT, where R is the reflux ratio, P0 is the pressure of the column, FD is the distillate flux, FS1−2 represent the side fluxes, FL1 is the liquid interconnection flow, FV2 is the vapor interconnection flow, NF is the feed stage number, NP0 is the first prefractioner tray location, NP is the last prefractioner tray location, NS1−2 represent the side stream tray locations, and NT is the number of column stages.Discrete and continuous variables:R,P0,FD,FS1,FS2,FL1,FV2,NF,NP0,NP,NS1,NS2, and NT.ym≥xm, where ym is the resulting purity vector and xm is the required purity vector.Cabrera-Ruiz et al. [221]Heat integrated distillation sequence cost optimizationGlobal optimizationMinimize TAC, TAC = minCsi,pi;si∈S,p∈P, where TAC is the total annual cost, C is the cost, S is the set containing all possible sequences, P is the set containing all possible column pressures, si is a distillation sequence, and p is the pressure of the column.Discrete and continuous variables:si and p.No information available about constraints.Lashkajani et al. [235]Water and energy minimization in the thermal cracking processGlobal optimizationMinimizeEC−minEC/maxEC−minEC + WU−minWU/maxWU−minWU,EC=EFuelgas+EFueloil/Product, WU=Reposition water/Product, where EC is the energy consumption, WU is the water use, EFuelgas is the gas energy consumption at the boiler, EFueloil is the oil energy consumption at the boiler, Product is the amount of produced product, Reposition water is the feed water for the cooling tower and steam grid, min is the minimum value, and max is the maximum value.Naphtha mass flow rate, ratio of steam dilution, temperature in the transfer line exchange, temperature of high-pressure steam, steam induction to turbines, distillate to feed ration column, temperature after transfer line exchange, and total naphtha mass flow rate.Process and utility constraints related to production and quality requirements, equipment performance, and resource characteristics.Boto et al. [239]Separation cascade design for multi-component mixturesGlobal optimizationMinimizec/∑Zn/ayk,P∗−yk,p+bxk,W∗−xk,W, where ∑Zn is the total flow rate of the cascade, yk,p∗ is the minimum concentration in the product, xk,W∗ is the maximum concentration in the waste, yk,p is the actual concentration in the product, xk,W is the actual concentration in the waste, and a,b,c are problem-related coefficients.Number of centrifuges and feed stage.Minimum concentration of componentk in the product yk,P≥yk,P∗. Maximum concentration of component k in the waste xk,W≥xk,W∗.Mansourzadeh et al. [244]Diesel blending problemGlobal optimizationMaximize∑j=1nQPj, where QPj is the final product flow rate and n is the total number of products.Component flow ratesQij.(1) Quantity conservation constraints:QFu=∑u′∈Uu∑s∈Su′,uQu′,u,s, QFu<QFumax, Qu,s=du,s.QFu,s∈Su, Qu,s=∑u′∈Us,uQu,s,u′, where QFu is the feed flow rate of unit u, Qu,s is the flow rate of stream s from unit u, Qu′,s,u is the flow rate of stream s coming from units u′ reaching unit u, Qu,s,u′ is the flow rate of stream s coming from unit u reaching unit u′, du,s is the yield level of stream s from unit u, Su is the set of streams from unit u, Su′,u is the set of streams leaving u′ units and reaching u unit, Uu is the set of units with u destination, and max is the maximum value. (2) Quality constitutive relation constraints: ∑i=1mai,j,kQi,j/∑i=1mQi,j=qj,k, j=1,2,…,n;k=1,2,…,l1, where ai,j,k is a component property, qj,k is a final product property, m is the number of component streams, j is the number of final products, and l1 is the number of linear properties. (3) Non-linear quality constitutive relation constraints: qj,l1+1=fQi,j,QPj,ai,j,l1+1, qj,l1+2=fQi,j,QPj,ai,j,l1+2, qj,l1+3=fQi,j,QPj,ai,j,l1+3.Gao et al. [242]Table 10
Problem description of life cycle analysis.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferenceLife cycle cost and CO2-eq emission optimizationMulti-objective optimizationMinimize LCC. MinimizeCO2−eq. emissions. LCC=IC+PVrplcmnt+PVOM&R+PVenergy. LCC is the life cycle cost. CO2−eq. represents the equivalent CO2 emissions. IC represents the initial investment costs. PVrplcmnt represents the present replacement cost value. PVOM&R is the present energy cost value. PVenergy is the present operation, maintenance, and repair cost value.Design parameters including construction materials of the roof, walls, windows, ceiling and floor, and the glazing type.Structure-related parameters and climate regions.UnconstrainedUnconstrainedFesanghary et al. [197]Asadi [230]Life cycle cost of the hybrid photovoltaic-diesel generator-battery systemGlobal optimizationMinimize LCC, LCC=∑m∈PV,DG,Battery,Conv/InvCCnpv,m+MCm+CDF, where LCC is the life cycle cost, CC is the annual capital cost, npv is the net present value, MCm is the annual maintenance cost, CDF is the annual cost of fuel consumption, PV is the set of photovoltaic panels, DG is the set of diesel generators, Battery is the set of batteries, and Conv/Inv is the set converter/inverters.Three decision variables: the number of batteries, the number of panels, and the diesel fuel consumption.NPV≥0, NBattery≥0, PPVt+PBt+PDGt≥Pdt, SOCmin≤SOCt≤SOCmax, SOCmin=1−DODSBattery, where NPV is the number of photovoltaic panels, NBattery is the number of batteries, PPV is the solar power, PB is the charge/discharge battery power, PDG is the diesel generation power, Pd is the total load demand, SOC is the battery bank charge, DOD is the maximum depth of the battery discharge, SBattery is the nominal battery capacity, t is the time, min is the minimum value, and max is the maximum value.Maleki et al. [236]Table 11
Problem description of water distribution system design.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferencePump-included water distribution system designGlobal optimizationMinimizeZ, Z=∑i=1NPipefPipeDi∗Li + ∑k=1NPumpfPump1QkRated,HkRated+CPV∗fPump2Qk,Hk,ηk, where Z is the total cost, Di is the pipe diameter, Li is the pipe length, fPipe is the pipe cost function, fPump1 is the pump capital cost, fPump2 is the annual pumping energy cost, QkRated is the rated pump discharge, HkRated is the rated pumping head, CPV is the present worth factor, Qk is the pump discharge, Hk is the pumping head, ηk is the pump efficiency, NPipe is the number of pipes, and NPump is the number of pumps.Pipe diameters and pump size.A penalty function was used to handle water pressure constraints:Hjmin≤Hj≤Hjmax, where Hj is the node pressure head, min is the minimum value, and max is the maximum value.Geem [216]Irrigation system design costGlobal optimizationMinimize Cost, Cost=∑i=1NCCDi+CMDi+CPDiLi, where Cost is the total cost, CC is the construction cost, CM is the maintenance cost, CP is the pipe material cost, Di is the pipe diameter, Li is the pipe length, and N is the number of pipes.Pipe diametersA penalty function was used to handle the pressure head and pipe velocity:hjmin≤hj≤hjmax, vimin≤vi≤vimax, where hj is the node pressure head, vi is the pipe velocity, min is the minimum value, and max is the maximum value.Yoo et al. [232]Pump scheduling problemMulti-objective optimizationMinimizeCT=∑t=124Ct∗Et, MinimizeCT=∑t=124NSWt, CT is the total energy cost, C is the unit cost, E is the consumed energy, t is the time, NSW is the number of pump switches.N number of active pumps at time t.TNSWp≤SW, Hmin≤Ht≤Hmax, 0≤Qpt≤Qpmaxt, Ht=24≥Ht=0, where TNSWp is the number pump switches, SW is the switch number per pump, H is the tank water level, Qp is the pump flow, min is the minimum value, and max is the maximum value.De Paola et al. [237]Water distribution system designMulti-objective optimizationMinimizeF1=∑in1.1∗Di1.24∗Li, MaximizeF2=αc=Pc/σc. F1 is the system cost design, F2 is the system robustness, Di is the pipe diameter, Li is the pipe length, n is the number of pipes, αc is the system robustness index, Pc is the mean stochastic nodal pressure at the critical node, σc is the standard deviation stochastic nodal pressure at the critical node.Pipe sizesHj≥17.6,j=1,…,m, where Hj is the pressure head and m is the number of consumer nodes.Jung et al. [238]Re-chlorination injection for water distribution networkGlobal optimizationMinimizef=∑i=1nMi, where f is the total disinfectant dose, Mi is the disinfectant mass, and n is the number of nodes.Booster injection mass rates for every node.(1) Residual chlorine requirements:Clmin≤Cli≤Clmax,i=1,…,n, where Cli is the concentration of residual chlorine in the consumer node, min is the minimum value, and max is the maximum value. (2) Re-chlorination facilities: N=1,…,n, where N is the number of re-chlorination facilities to install. (3) Hydraulic constraints: ∑Qin−∑Qout=Qe, where Qin is the incoming node flow, Qout is the outgoing node flow, and Qe is the external node demand. (4) Energy conservation: ∑k∈LooplΔHk=0,∀l∈NL, where ΔHk is the pipe head loss and NL is the total number of loops.Yoo et al. [240]Table 12
Problem description of other applications in chemical engineering.
ApplicationType of optimization problemObjective functionOptimization variablesConstraintsReferenceCombined heat and power economic dispatch problemGlobal optimizationMinimizeC, C=∑i=1npcipi+∑j=1nccjhjpj+∑k=1nhckhk, where C is the power and heat total production cost, ci is the unit production cost, pi is the unit power generation, hj is the unit heat production, np is the number of power units, nc is the number of cogeneration units, and nh is the number of heat units.pi,pj,hj, and hk∑i=1nppi+∑j=1ncpj=PD, ∑j=1nchj+∑k=1nhhk=HD, where PD is the system power demand and HD is the system heat demand.Vasebi et al. [214]Chemical equation balancingGlobal optimizationMinimize∑j=1NHjxi, Hjxi=hjxi4/∑i=1Mhjxi4, where Hj is the jth normalized constraint associated with xi, xi is the stoichiometric coefficient of the corresponding reactant or product, hjxi is the jth constraint associated with xi, N is the number of elements, and M is the number of different reactants and products.Stoichiometric coefficients of reactants and productsxi.hjxi=∑i=1Majixi,j=1,…,N, where aji is the number of atoms of the element j in each reactant or product i.Zou et al. [219]Pollution source identification for groundwaterGlobal optimizationMinimize∑k=1nk∑i=1nobwcikCestik−Cobsik2, wcik=1/Cobsik+η2, where wcik is the corresponding weight, Cestik is the estimated concentration, Cobsik is the observed concentration, nob is the number of observation locations, nk is the number of observation time periods, and η is a non-domination constant.Fluxes and locations of pollution sources.C=fX,Y,q, where C is a vector containing estimated concentrations, X,Y are vectors containing source locations, q is a vector containing source fluxes, Xmin≤X≤Xmax, Ymin≤Y≤Ymax, qmin≤q≤qmax, min is the minimum value, and max is the maximum value.Jiang et al. [227]Kinetic parameter estimation problemGlobal optimizationMinimize1/m∑i=0m−1|xLi−x^Li|/x^Li, where m is the number of samples, xL is the condensation product fraction, x^L is the fraction obtained from x^L=KLP0e−ELP/T/nL1−1−xnL + KWP0KWLP0e−EWP+EWLP/T/nW−KWLP0e−EWLP/T ×1/KWLP0e−EWLP/T1−1−xKWLP0e−EWLP/T+1/nW1−xnW−1, x is the mole fraction, T is the temperature, and KLP0,ELP,KWP0,KWLP0,EWP,EWLP,nW, and nL are the adjustable parameters of the model.Model parametersKLP0,ELP,KWP0,KWLP0,EWP,EWLP,nW, and nL.No information available about constraints.Ma et al. [228]Inverse problem of chemical kineticsGlobal optimizationMinimize∑i=1ndat∑j=1Nxijcalc−xijexp, where xij is the mole fraction, calc is the calculated value, exp is the experimental value, ndat is the number of measured points, and N is the number of substances in the reaction.Activation energies, pre-exponential multipliers, constant parameter, and denominator order.UnconstrainedEnikeeva et al. [251]Multiple effect evaporator designGlobal optimizationMinimize WRF, WRF=∫Wt×Rtdt, Rt=ft,xit,xi′t,Tit,Ti′t, where WRF is the weighted residual function, W is the weight function, R is the resulting residual function from the composition and temperature approximations and their derivatives, xi is the exit liquor concentration, xi′ is the derivative of the exit liquor concentration approximation function, Ti is the temperature of the generated vapor, Ti′ is the derivate of the generated vapor temperature approximation function, and t is the time.Fourier series coefficientsa0,am,bm,c0,cm, and dm.A penalty function method was used to handle constraints:Ti≥Ti+1i=1,…,6,T7≥52T8.Yadav et al. [250]


### 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 compositionz, 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 = 10∗nvar, where nvar 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 = UBi − LBi, PAR = 0.75, and HMCR = 0.5 for HS and GHS and bwmin = 0.001, bwmax = UBi − LBi, PARmin = 0.5, PARmax = 0.95, and HMCR = 0.5 for IHS. For all cases, HMS was defined as 10∗nvar. Two termination criteria were used: a maximum number of iterations without improvement and itermax. 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 itermax = 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. bwi = UBi − LBi, PAR = 0.75, HMCR = 0.25, and HMS = 10∗nvar 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 = 10∗nvar, HMCR = 0.5, PAR = 0.5, and itermax = 20000∗nvar. 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-CO2 + α-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, bwini = 0.5, bwmin = 1000, bwmax = 2, bwsat = 1000, and amplitude constant Cbw = 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: itermax = 5E104, 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 = 10∗nvar. 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: bwi minimum and maximum values were defined in terms of the optimization variables bounds, PARmin = 0.1, PARmax = 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 bwi 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 itermax). 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 (CO2-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, CO2 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 CO2 emissions with an increase in LCC. Asadi [230] also minimized the LCC and CO2-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 itermax = 1,000, bwmin = 0.01, bwmax = 1, PARmin = 0.1, PARmax = 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 1012 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 m3 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 m3 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 bwmin = 1/(30 (UBi − LBi)), bwmax = 1/(20 (UBi − LBi)), PARmin = 0.01, PARmax = 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: itermax = 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.


## 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 compositionz, 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 = 10∗nvar, where nvar 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 = UBi − LBi, PAR = 0.75, and HMCR = 0.5 for HS and GHS and bwmin = 0.001, bwmax = UBi − LBi, PARmin = 0.5, PARmax = 0.95, and HMCR = 0.5 for IHS. For all cases, HMS was defined as 10∗nvar. Two termination criteria were used: a maximum number of iterations without improvement and itermax. 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 itermax = 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. bwi = UBi − LBi, PAR = 0.75, HMCR = 0.25, and HMS = 10∗nvar 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 = 10∗nvar, HMCR = 0.5, PAR = 0.5, and itermax = 20000∗nvar. 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-CO2 + α-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, bwini = 0.5, bwmin = 1000, bwmax = 2, bwsat = 1000, and amplitude constant Cbw = 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: itermax = 5E104, 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 = 10∗nvar. Global performance results placed HS in the third position of tested stochastic optimizers.


## 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 compositionz, 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 = 10∗nvar, where nvar 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 = UBi − LBi, PAR = 0.75, and HMCR = 0.5 for HS and GHS and bwmin = 0.001, bwmax = UBi − LBi, PARmin = 0.5, PARmax = 0.95, and HMCR = 0.5 for IHS. For all cases, HMS was defined as 10∗nvar. Two termination criteria were used: a maximum number of iterations without improvement and itermax. 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 itermax = 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. bwi = UBi − LBi, PAR = 0.75, HMCR = 0.25, and HMS = 10∗nvar 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 = 10∗nvar, HMCR = 0.5, PAR = 0.5, and itermax = 20000∗nvar. 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-CO2 + α-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, bwini = 0.5, bwmin = 1000, bwmax = 2, bwsat = 1000, and amplitude constant Cbw = 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: itermax = 5E104, 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 = 10∗nvar. 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: bwi minimum and maximum values were defined in terms of the optimization variables bounds, PARmin = 0.1, PARmax = 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 bwi 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 itermax). 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 (CO2-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, CO2 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 CO2 emissions with an increase in LCC. Asadi [230] also minimized the LCC and CO2-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 itermax = 1,000, bwmin = 0.01, bwmax = 1, PARmin = 0.1, PARmax = 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 1012 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 m3 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 m3 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 bwmin = 1/(30 (UBi − LBi)), bwmax = 1/(20 (UBi − LBi)), PARmin = 0.01, PARmax = 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: itermax = 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.



---

*Source: 1928343-2022-10-07.xml*