A Tandem Running Strategy-Based Heat Transfer Search Algorithm and Its Application to Chemical Constrained Process Optimization

Constrained optimization problems (COPs) are widely encountered in chemical engineering processes, and are normally defined by complex objective functions with a large number of constraints. Classical optimization methods often fail to solve such problems. In this paper, to solve COPs efficiently, a two-phase search method based on a heat transfer search (HTS) algorithm and a tandem running (TR) strategy is proposed. The main framework of the MHTS–TR method aims to alternate between a feasible search phase that only examines feasible solutions, using the HTS algorithm, and an infeasible search phase where the treatment of infeasible solutions is relaxed in a controlled manner, using the TR strategy. These two phases play different roles in the search process; the former ensures an intensified optimum in a relevant feasible region, whereas the latter is used to introduce more diversity into the former. Thus, the ensemble of these two complementary phases can provide an effective method to solve a wide variety of COPs. The proposed variant was investigated over 24 well-known constrained benchmark functions, and then compared with various well-established metaheuristic approaches. Furthermore, it was applied for solving a chemical COP. The promising results demonstrate that the MHTS–TR approach is applicable for solving real-world COPs.


Introduction
Nowadays, many real-world chemical engineering processes are defined by complex objective functions with a large number of constraints [1]. The optimization problems that contain several constraints are described as constrained optimization problems (COPs) [2]. These problems are generally characterized by their different types, such as linear, nonlinear, polynomial, quadratic, cubic, etc. Due to the complexity of highly constrained chemical processes, new generation optimization methods need to be found, as classical methods often fail to solve COPs efficiently. Thus, various metaheuristic algorithms (MHAs) have been developed, modified, and applied significantly to optimize a wide variety of COPs [2][3][4].
Various approaches have been applied as constraint-handling techniques to deal with COPs during the search course, as reported in the surveys [5,6]. One of the most popular approaches are penalty-based methods [7][8][9], which can be simply classified into static and self-adaptive methods. The static methods deal with the infeasible solutions, by transforming a COP into an unconstrained problem. However, defining the penalty parameter values is not straightforward. Self-adaptive penalty methods modify the penalty term value throughout the search course, such as the adaptive penalty method (APM) [10], an effective penalty-based method that automatically calibrates the infeasible surface throughout evolution. However, it may lose feasible solutions during the search course. A set of three feasibility rules called Deb's rules [11] is a popular selecting candidate solutions technique, using pre-established priority rules which separately rank the solutions depending on their constraint violation values and objective function values, respectively. These rules state that: (i) between two feasible solutions, the solution that has a superior value of the objective function is preferred; (ii) between feasible and infeasible solutions, the feasible one is preferred; (iii) between two infeasible solutions, the one with smaller constraint violation values is preferred. These rules can be incorporated with several MHAs, since they do not often require any additional parameters [12]. A variant of the feasibility rules for diversity maintenance was proposed by Mezura and Coello [13]. It is one of the simplest mechanisms that allows a set of infeasible solutions to remain in the search. In this variant, these solutions are adjacent to the feasible domain; moreover, they have a great objective function value and the lowest sum of constraint violation. These solutions are selected from the offspring or the parent's population with a 50% probability.
When solving complex COPs, simply finding a feasible solution is not a straightforward task. Several types of research on highly constrained optimization problems [14][15][16] have reported that the consideration of infeasible solutions during the search course, instead of limiting the search process to feasible regions, may help to better explore the search space. As it can be noted from the above-mentioned works of literature, the main goal of their work is the focus on the addition of evolutionary mechanisms within the feasible and infeasible regions during the search process. The idea of preserving the infeasible solutions close to the feasible region allows for finding an optimum on the boundary of the feasible region of the search space [17]. The optimal solution of COPs is often located in the interior or near the boundary of the feasible region. Hence, the individuals that are far away from the feasible region are practically unhelpful to the optimization of the population, whereas the individuals that are close to the feasible region may contain favorable information that can help the population in searching for an optimum on or near the boundary of the feasible region. Therefore, based on the above observation, it is important to consider how to properly preserve the infeasible individuals close to the feasible region during the search course, for the sake of guiding the search of the feasible region toward the global optima of COPs.
The heat transfer search (HTS) algorithm [18] is a novel population-based method inspired by the natural laws of thermodynamics and heat transfer. Its basic framework mainly includes three main phases such as conduction, convection, and radiation. This algorithm belongs to the family of new generation MHAs, and is considered one of the most competent optimization approaches in comparison with other MHAs, as shown in the literature [18]. Although it is a relatively new method, a number of its variants have been developed, modified, and applied significantly to solve a wide variety of optimization problems [19][20][21][22][23][24][25][26]. However, the HTS-based variants have not been applied for handling chemical COPs. Therefore, we aim to extend the actual applications of this MHA to handle such problems.
The major contributions of our paper are described below: (1) A novel approach with two search phases called MHTS-TR is proposed through integrating the multiple HTS algorithm and the TR strategy. The ensemble of these two complementary phases can provide an effective algorithm for solving various COPs; (2) The effectiveness of the new variant is tested through a set of 24 constrained benchmark problems, and the simulation results are compared with those of other competitors; (3) The MHTS-TR approach is used to handle a real-world chemical COPs. Additionally, the simulation results obtained on this problem are compared with those of different approaches existing in the literature. To the best of our knowledge, this paper presents the first attempt for applying an HTS-based method to handle a chemical COP.
The rest of our work is organized as follows: the main formula of COPs is defined in Section 2; the main theoretical principles of the TR strategy and HTS method are described in Section 3; the new variant is explained in detail in Section 4; in Section 5, certain experimental investigations and comparisons are conducted, and the proposed MHTS-TR approach is used for solving a real-world COP in Section 6. In Section 7, our final conclusions are summarized.

Problem Statement
In general, the mathematical model of a COP can be described as follows, where the main goal is to optimize the objective function, represented as f (x): (1) where f (x) represents the fitness function; x ∈ Ω ⊆ S indicates the n-dimensional solution vector, x i denotes the ith dimensional component of x; S ∈ R n indicates the solution space determined by the upper and lower bounds (x max = [x 1 max , . . . , x i max , . . . , x n max ] and x min = x 1 min , . . . , x i min , . . . , x n min ) of the solution vector x; Ω represents the feasible region of dimension n; g j (x) ≤ 0 indicates the inequality constraint; h j (x) = 0 denotes the equality constraint, and l and p are defined as the number of inequality and equality constraints, respectively. Due to the constraints shown in Equation (1), two disjoint subsets (feasible and infeasible) constitute the search domain. The feasible domain is defined by the regions where all p constraint functions of equalities and inequalities are satisfied. Thus, the solutions x belonging to the feasible region and infeasible region are classified as feasible and infeasible candidate solutions, respectively.
In general, the constraint-handling techniques can be classified either as indirect, when both feasible and infeasible candidate solutions are considered along the search, or as direct, when only the feasible candidate solutions are employed. The penalty method is the most common indirect approach applied with MHAs to penalize the infeasible solutions. In this method, when x is an infeasible solution, its objective function is penalized by adding a penalty term, which depends on the constraint violation. When solving COPs, in addition to calculating the fitness value of the population according to the function f (x), it is also necessary to evaluate the constraint violation. Generally, the violation degree of a member x to the jth constraint can be expressed as follows: Here, the absolute value of the equality constraint function ( h j (x) ) can be treated as an inequality given by G j (x) = max 0, h j (x) − δ , where δ is a small positive value.
The general form of the penalty function (p (x)) and the corresponding evaluation function (eval (x)) can be described as follows [1]: where β and C are generalized dynamic or static coefficients, chosen according to the applied technique; F and U represent the feasible and infeasible spaces, respectively. When handling COPs, p(x) is usually used to evaluate the infeasibility of the population.

Heat Transfer Search (HTS) Algorithm
The HTS algorithm is a relatively new population-based approach that belongs to the family of MHAs. It is inspired by the natural laws of thermodynamics and heat transfer; [18] declares that "any system usually attempts to attain an equilibrium state with the surroundings" [18]. It has been reported that the HTS algorithm mimics the thermal equilibrium behavior of the systems by considering three heat transfer phases, including the conduction phase, convection phase, and radiation phase [18]; each phase plays a crucial role in establishing the thermal equilibrium and attaining an equilibrium state. Similarly to other MHAs, this algorithm starts with a randomly initialized population, and the population is considered as a cluster of the system's molecules. These molecules aim to attain an equilibrium state with the surroundings through the three phases of heat transfer, by interacting with each other and with their surrounding environment. In the basic HTS algorithm, the population members are only updated through one phase of the three heat transfer phases in each iteration. The selection process for which of the three phases to be activated for updating the solutions in the particular iteration is carried out by a uniformly distributed random number R. This random number is generated in the range [0, 1], randomly, in each iteration to determine the phase that must be selected. In other words, the population members undergo the conduction phase when the random number R varies between 0 and 0.3333, the radiation phase when the random number R varies between 0.3333 and 0.6666, and the convection phase when the random number R varies between 0.6666 and 1. The greedy selection technique is the main selection mechanism for newly generated solutions in the HTS algorithm; this technique states that only new updated solutions which have a superior objective value will be accepted, and the solutions with an inferior objective value will be subsequently substituted by the best solutions. Hence, by comparing the difference between the current solution and the elite solutions, the greatest solution can be finally achieved. In the basic HTS algorithm, the main search process is performed by the elementary operations of the three heat transfer phases (conduction, convection, and radiation); the basic principle of each phase is briefly described in the subsequent subsections. The overall flow-chart of the original HTS method is illustrated in Figure 1.

Conduction Phase
The fundamental principle of the conduction phase in the original HTS algorithm is concerned with the influence of conduction among the substance's molecules, and the outcome is that the heat transfer takes place. Thus, by aiming to attain a state of thermal equilibrium through the conduction influence, the molecules that have higher energy transfer the heat to the molecules that have lower energy. In other words, the system here aims to neutralize the thermal imbalance among the substance's molecules. In the conduction phase of the HTS algorithm, a potential new solution is updated by utilizing the following formulas: where maxFEs denotes the maximum function evaluations; CDF indicates the conduction factor which is set to a value of 2 for balancing the exploration and exploitation in this phase; FE indicates the function evaluations; X new j,i denotes the newly updated solution; j = 1, 2, . . . , n; k indicates a randomly chosen solution; j = k, where k ∈ (1, 2, . . . , n); i represents a randomly chosen decision variable, and i ∈ (1, 2, . . . , m); R indicates the random variable that varies in the interval [0, 0.3333], and r i indicates a random number that varies in the interval [0, 1]. Here, r i and R 2 are in correlation with the conductance elements of the Fourier's law of heat conduction [18].

Convection Phase
The fundamental principle of the convection phase in the basic HTS algorithm is that due to the convection effect between the surrounding temperature and the system, heat transfer occurs. Thus, by attaining a state of thermal equilibrium through the convection effect, the mean temperature of the system interacts with the adjacent surrounding temperature to neutralize the thermal imbalance. In the convection phase, the mean temperature of the system is represented by the mean of the population members, which is marked as X ms , whereas the surrounding temperature is represented by the best solution, which is marked as X s . Therefore, a new solution is updated in this part of the algorithm using the following equation: where COF represents the convection factor, and is set to a value of 10 for balancing the exploration and exploitation in this phase; FE denotes the function evaluations; TCF indicates the temperature change factor; X new j,i denotes the new updated solution; j = 1, 2, . . . , n; i = 1, 2, . . . , m; r i is a random number which varies in the interval [0, 1], and R is the random variable that varies in the interval [0.6666, 1]. Here, r i and R are in correlation with the convection elements of Newton's law of cooling [18]. The value of TCF varies randomly in the range [0, 1] in the first part of this phase, whereas it varies either as 1 or 2 in the second part of this phase.

Radiation Phase
The fundamental principle of the radiation phase in the original HTS algorithm is that the system attempts to neutralize the thermal imbalance by interacting within the system itself (for example, the other solution), or with the surrounding temperature (for example, the best solution). Here, the system aims to establish a state of thermal balance. In the radiation phase of the HTS algorithm, a potential new solution is updated by utilizing the subsequent formulas: where RDF denotes the radiation factor whose value is set at 2 for balancing the exploitation and exploration in this phase; FE indicates the function evaluations; X new j,i denotes the new updated solution; j = 1, 2, . . . , n; i = 1, 2, . . . , m; k represents the randomly selected solutions; j = k, and k ∈ (1, 2, . . . , n); r i indicates a random number which varies in the interval [0, 1], and R denotes the random variable which varies in the interval [0.3333, 0.6666]. Here, r i and R are in correlation with the radiation parameters of the Stefan-Boltzmann law [18].

Proposed MHTS-TR Algorithm
In this section, by integrating the tandem running (TR) strategy into the basic HTS algorithm, a novel variant of the HTS algorithm named MHTS-TR is presented. The proposed variant is described in detail in the following subsections.

Tandem-Running Strategy
In nature, animals acquire information socially, and comprehensively utilize such information to allow them to exist in large, collaborative societies. The behavior of energetically sharing personal information with other members of a social group is referred to as cooperative communication [27]. Tandem running (TR) is a form of recruitment seen in some groups of ants, in which one member guides another to a crucial resource, such as a superior nest location. For instance, the Temnothorax Albipennis ant species use tandem-run recruitment, where one ant acts as a leader and guides a single follower from the nest to food; thus, the follower finds food more rapidly than when seeking alone. The tandem leader is familiar with the location of the food resource whereas the follower is naive. The course and speed of the leader are governed by the follower, who taps the leader's legs and gaster with its antennae [28]. The following ant acquires knowledge regarding the placemarks and surroundings en route, and thus learns to take a more linear route, usually returning to its nest quicker than before [28]. Generally speaking, when ant colonies are suddenly placed in a complex and unfamiliar environment, this information sharing behavior commonly appears. Thus, this behavior can also be mimicked in complex environments for exploring useful formation.

Outline of the MHTS-TR Algorithm
Inspired by the cooperative exploration behavior of the Temnothorax Albipennis ant species in unknown areas, described as a "move-to-improve" method, this strategy is incorporated into the HTS algorithm, and a new variant named MHTS-TR has been proposed. The overall scheme of the MHTS-TR algorithm is shown in Figure 2. For this approach, we divided the population members into two categories (feasible and infeasible) at each iteration, depending on constraint violation. The members that did not violate the constraints were termed as "leaders", and were evolved within the feasible region by the HTS algorithm. By contrast, the violated members were termed as "followers", and were further classified into two parts (X HV and X SV ) depending on their violation degree. X HV represented the member with a higher degree of violation, with a position often further away from the feasible region. Thus, a member was randomly selected from the feasible region to be the leader, and the followers then moved toward the neighborhood of this leader to search within the feasible region. As a result, the members within the infeasible region that did not contribute to the population were moved toward the feasible region. Meanwhile, as each follower randomly selected its leader, the population density within the feasible region increased evenly and, hence, increased the diversity of the population within the feasible region. On another hand, X SV represented the member with a relatively lower degree of violation, which was considered as it was almost close to the feasible region. It selected the nearest member that was located in the feasible region to be its leader, and moved towards it; hence, the boundaries of the feasible region were gradually searched by approximating towards the leader. In this way, the members with infeasible information nearby the boundaries were utilized to explore the superior areas that were hidden nearby the boundaries of the feasible region. Therefore, due to the methods used by X HV and X SV to select their respective leaders being carried out through random selection and distance judgment, which was irrespective of the fitness value, there was no issue of the members being overly concentrated around the global optimal member. Thus, the non-connectivity in the feasible region did not affect the distribution of the population in each feasible region. By alternating between these two complementary phases, the MHTS-TR method was expected to explore various zones of the search space without being easily trapped in a local optimum. The moving strategies of X HV and X SV are shown in Figure 3.

The Overall Process of MHTS-TR Method
Firstly, we assumed that the population M was the number of members that searched in an n-dimensional space (S ∈ R n ), and Ω ⊆ S was the feasible region of the solution space. At the beginning of the kth iteration, the distance matrix (Dis k = Dis 1 k , . . . , Dis i k , . . . , Dis M k ) for all the members was calculated, in which Dis i k was an m-dimensional vector that represented the distance between the member i and other members, and where dis ij k was the Euclidean distance between the member i and member j (1 ≤ j ≤ M and j = i).
Secondly, all members of the population were evaluated as to whether they were within the feasible region or not, and then divided into two sub-populations, accordingly. We assumed that Fb was the number of members within the infeasible region, where 0 ≤ Fb k ≤ M, and the number of members within the feasible region was M − Fb k . These members were considered as "leaders", and would evolve according to the basic rules of the conduction, convection, and radiation phases of the HTS algorithm.
The degree of constraint violations for all the members in the infeasible region was calculated G(x) = [G(x 1 ), . . . , G(x i ), . . . , G(x Fb )], i = 1, 2, . . . , Fb, and then G(x) was sorted in descending order. A number of members with the highest degree of violation (denoted by SN) were selected to be X HV , and the other members were selected to be X SV . Here, the SN was jointly determined by the number of members within the infeasible region (Fb) and the X HV ratio (ps k ) at the kth iteration. It can be calculated as follows: where the operational sign " Processes 2021, 9, x FOR PEER REVIEW 9 of 23

The overall Process of MHTS-TR Method
Firstly, we assumed that the population M was the number of members that searched in an n-dimensional space (  n SR ), and S was the feasible region of the solution space. At the beginning of the kth iteration, the distance matrix (

Dis Dis Dis Dis
) for all the members was calculated, in which k i Dis was an m-dimensional vector that represented the distance between the member i and other members, and 1 , , ,where k ij dis was the Euclidean distance between the member i and member j (1 ≤ j ≤ M and j ≠ i).
Secondly, all members of the population were evaluated as to whether they were within the feasible region or not, and then divided into two sub-populations, accordingly. We assumed that Fb was the number of members within the infeasible region, where 0  k Fb M , and the number of members within the feasible region was − k M Fb . These members were considered as ''leaders'', and would evolve according to the basic rules of the conduction, convection, and radiation phases of the HTS algorithm. The degree of constraint violations for all the members in the infeasible region was where the operational sign "   " means round number; increased with the increase in iterations throughout the whole process of calculation, meaning that the MHTS-TR algorithm had a high tolerance level for the violated members at the beginning of the calculation. In other words, most of the members within the infeasible region were moved towards the feasible region using the approximation manner, whereas in the later stage of the calculation, the tolerance of violated members was gradually tightened by the increasing the ratio of XHV. Most of them were directly updated into the interior of the feasible region for searching. The detailed evolutionary mechanisms of the XHV and XSV are given as follows: (1) The moving strategy of XHV: " means round number; ps min and ps max are the initial and terminal values of ps, respectively, where 0 ≤ ps min < ps max ≤ 1. The value of ps k increased with the increase in iterations throughout the whole process of calculation, meaning that the MHTS-TR algorithm had a high tolerance level for the violated members at the beginning of the calculation. In other words, most of the members within the infeasible region were moved towards the feasible region using the approximation manner, whereas in the later stage of the calculation, the tolerance of violated members was gradually tightened by the increasing the ratio of X HV . Most of them were directly updated into the interior of the feasible region for searching. The detailed evolutionary mechanisms of the X HV and X SV are given as follows: (1) The moving strategy of X HV : We assume that, a member i is selected to be X HV at the kth iteration, when its position is x k i = (x k i1 , x k i2 , · · ·x k in ) / ∈ Ω, then a member j will be randomly selected from the feasible region to be the leader of this member, so that the member i will be updated to a new position around the member j, where the position of the member j is x k j = (x k j1 , x k j2 , · · ·x k in ) ∈ Ω. Firstly, the leader (member j) will select the closest member g which is also within the feasible region according to the distance Dis j k , when its position is x k g = (x k g1 , x k g2 , · · ·x k gn ) ∈ Ω, and then the Euclidian distance between the member j and member g on each dimension is calculated NDis j k = [edis 1 , . . . , edis h , . . . edis n ], where edis h is the distance between the member j and member g on the h dimension, and its calculation formula is as below: Then, the X HV (member i) will update its position using the subsequent formula.
x k+1 where rand(n) is an n-dimensional random vector, it is uniformly distributed between 0 and 1, and the operator "×" means calculating the element-wise product of the two vectors.
(2) The moving strategy of X SV : We assume that, at the kth iteration, when a member i is selected to be X SV , it selects the nearest member j within the feasible region depending on the distance Dis i k between this member and other members, and then approximates towards the member j. Thus, it will update its position as follows: where c is the velocity factor of X SV , it is used to adjust the velocity of the X SV to approximate toward the feasible region. The overall process of this approach is illustrated by the flowchart diagram in Figure 4.

Numerical Experiments and Discussion
In this section, we outline how the overall effectiveness of the MHTS-TR approach was verified by a set of 24 well-defined COPs of Congress on Evolutionary Computation 2006 (CEC 2006) [29][30][31]. Moreover, comparisons between the new variant and several other well-established MHAs, such as differential evolution (DE), particle swarm optimization (PSO), biogeography-based optimization (BBO), artificial bee colony (ABC), teaching-learning-based optimization (TLBO), and the original heat transfer search (HTS) algorithm were conducted. These comparative approaches were examined against the considered benchmark problems noted previously in the literature [18]. Therefore, they were employed for comparison with the proposed variant; this is notable since a common experimental platform is required to make fair comparisons against the competitor algorithms. Hence, the population size (NP) was set at 50, and the maximum number of function evaluations (maxIter) was set to 240,000. Furthermore, the computational results obtained from 100 independent runs, such as the best value (Best), mean value (Mean), worst value (Worst), standard deviations (Std), and success rate (SR) were listed. To further quantify the performances of the competitors statistically, a Friedman rank test was also conducted on the mean and best solutions obtained by the comparative approaches. The constraint-handling technique used with the considered competitors was the static penalty method.

Benchmark Functions
The considered test functions offered challenges to the MHAs and can be classified into different types: non-linear, linear, quadratic, polynomial, and cubic functions. The detailed characteristics of each test function are illustrated in Table 1, and the detailed mathematical formulation of each test function with its constraints can be found in the literature [18].

Simulation Results and Analysis
The comparative results of each approach on C01 to C13 are displayed in Table 2. It can be observed from the results table that the MHTS-TR, HTS, and ABC algorithms markedly outperformed the other four algorithms on the results of the quadratic function (C01). Moreover, they produced a 100% success rate compared to the other algorithms on this function. In terms of the non-linear function (C02), although the ABC method produced superior results to the other approaches, all comparative results were considered, as they failed to attain the global optimum. For the polynomial test function (C03), the ABC approach produced superior results than other methods; however, the MHTS-TR and HTS methods obtained competitive performances, and achieved the global optimum on this test function. On the quadratic function (C04), all approaches produced equally great results except for the BBO algorithm. On the cubic function (C05), the MHTS-TR, HTS, and TLBO methods showed competitive results over the other methods, and the mean value acquired by the MHTS-TR was superior to that of the HTS and TLBO, respectively. On the cubic function (C06), the MHTS-TR, ABC, TLBO, HTS, and PSO approaches produced equally great values. On the quadratic function (C07), the mean value obtained by the MHTS-TR, HTS, ABC, and DE algorithms was better than those of the other competitive algorithms. On the non-linear function (C08) and polynomial function (C09), all the competitive methods (except the BBO) produced equally good results. On the linear function (C10), the mean value acquired by the MHTS-TR was superior to those of the other competitive algorithms. All the competitive algorithms (except the DE and BBO) produced equally great results on the quadratic function (C11). On the quadratic function (C12), all approaches (except the PSO) produced equally good solutions. All the competitive approaches failed to attain the global optimum value on the non-linear function (C13); however, the mean value obtained by the MHTS-TR and PSO approaches was better than those of other competitive algorithms.   Table 3 shows the Friedman rank test, which was conducted on the mean and best solutions obtained by the comparative methods on the functions C01 to C13. It can be seen from the result table that the MHTS-TR algorithm ranked first to obtain the mean solution, followed by the original HTS, ABC, and the rest of the other approaches, whereas it secured the second rank after the PSO algorithm to obtain the best solution on these test functions.
The results of the comparative algorithms on the functions C14 to C24 are shown in Table 4. It can be observed from the results table that all the methods failed to attain the global optimum value on the non-linear function (C14), which was a difficult function. However, the MHTS-TR algorithm produced a superior mean value compared to the other approaches. The mean result obtained by the MHTS-TR and DE algorithms outperformed the other algorithms on the quadratic function (C15). On the non-linear function (C16), it can be seen that all approaches (except the BBO) showed identical results and have succeeded to attain the global optimum value. On the non-linear functions (C17 and C19), and the quadratic function (C18), the mean results obtained by the TLBO, DE, and ABC algorithms were superior to the results of the other comparative algorithms. However, the MHTS-TR algorithm generated better results than the original HTS on these functions. On the rest of the linear functions (C20 to C24), due to C20, C22 and C23 having test problems as difficult functions, it can be seen that all the approaches failed to attain the global optimum value. However, the MHTS-TR method generated superior mean results compared to the other comparative approaches on the functions C20 and C23. The mean result and success rate obtained by the MHTS-TR algorithm on the function C21 was better than those of the other algorithms. On the function C24, it can be observed that all algorithms attained the global optimum value, and also shown identical results. Table 3. Friedman rank of C01-C13 for the solutions obtained by the MHTS-TR and other MHAs. (The experimental results of the competitors are derived from the literature [18]).

The Mean Solution
The   Table 5 shows the Friedman rank test, which was conducted on the mean and best solutions acquired by the comparative methods on the functions C14 to C24. It can be seen from the results table that the MHTS-TR algorithm ranked first to obtain the mean and best solutions, followed by the original HTS, TLBO, DE, and the other algorithms. Due to the appropriate trade-off between exploitation and exploration, our method was carried out by integrating different search mechanisms. The proposed method could alleviate the probability of local optima stagnation, and finally converge to a global feasible solution. Hence, it achieved the global optimum solution for more test problems compared to the original HTS method and the rest of the other competitive methods.

Further Discussion
To further analyze the advantages of the MHTS-TR algorithm in handling COPs, a comparison between the basic HTS and MHTS-TR on the test problem C10 was conducted. This comparison was carried out by counting the number of feasible solutions during the optimization course. To facilitate the comparison, the population size of both approaches was set at 100, and the maximum number of iterations was set to 1000. Figure 5 shows the number of feasible solutions obtained by both approaches in 1000 iterations. It can be seen that both algorithms found feasible solutions in the 10th iteration, and the number of feasible solutions was increased to the maximum value around the 100th iteration. In the later stage of calculation, the number of feasible solutions obtained by the HTS algorithm accounted for about 50% of the total number of members, and keptdecreasing at the later stages of the calculation. By contrast, the number of feasible solutions obtained by the MHTS-TR algorithm accounted for almost 70% of the total number of members, and was relatively stable. This implies that the improvements on the proposed MHTS-TR approach could greatly increase the population density and diversity, and improve the exploration and exploitation capabilities of the population in the feasible region.

Application to Chemical Constrained Process Optimization
In this section, we further assessed the overall performance of the MHTS-TR approach by applying it to a real-world chemical COP, which was the simplified alkylation process. Moreover, we compared the experimental results obtained by the MHTS-TR approach for this problem with those of other optimization approaches such as the α-based branch and bound method (αBB) [32], the constrained ant colony system (CACS) [33], cultural algorithm with evolutionary programming (CAEP) [34], branch and reduce optimization navigator (BARON) [35] and the basic HTS algorithm. Here, BARON version 21.1.13 with default options was used. The first three methods were tested on this chemical COP previously in the literature; hence, they were employed for comparison with the new variant. To be consistent with the considered competitors, the maximum number of iterations (maxIter) was set to 1000, and the population size (NP) was set at 50. A detailed description of this chemical process is given in the subsequent sub-section.

Problem Description
The simplified alkylation process was firstly discussed in detail by Bracken and McCormick [36]. The flowsheet of this process is shown in Figure 6. In this process, an olefin feed of 100% butane, a 100% isobutane make-up stream, and a pure isobutane recycle were infused into a reactor. Then, these components reacted with an acid catalyst, and some spent acid was discharged out of the reactor. Moreover, the product in the reactor was transported into a fractionator, and separated into isobutane and alkylate products. Bracken and McCormick [36] proposed an optimizing model for the alkylation process, and the objective was to maximize the profit. Dembo [37] transformed this process into a model with 7 variables and, subsequently, a slightly modified version of this model was proposed by other researchers in the literature [32]. The variables of this process are listed in Table 6. Motor octane number - External isobutane-to-olefin ratio % x 7 F−4 performance number - The profit-maximization problem for this process takes the form as follows [32]: g 12 = 1020.4082x 4 x 2 + 1.2244898x 3 x 4 − 100, 000x 2 ≤ 0, g 13 = 6.25x 1 x 6 + 6.25x 1 − 7.625x 3 − 100, 000 ≤ 0,

Simulation Results and Discussion
The best values of the profit objective obtained by each approach are displayed in Table 7. The mean execution times of the original HTS and the MHTS-TR approaches on this problem were as follows: the original HTS algorithm spent 12.38 s, and the MHTS-TR approach spent 14.21 s to get these results. It can be observed from Table 7 that the optimization result obtained by the MHTS-TR approach was better than those of the basic HTS, BARON, CACS (δ = 5 × 10 −6 ), and CACS (δ = 0) methods. Moreover, the results obtained by the αBB, CAEP, and CACS (δ = 5 × 10 −4 ) methods were superior to that of the MHTS-TR, but the MHTS-TR method did not violate any constraint. The violations of constraints for the best solutions generated by all algorithms are listed in Table 8. It can be observed from the results table that the αBB, CAEP, CACS (δ = 0), and CACS (δ = 5 × 10 −4 ) approaches violated at least one of the constraints, whereas the CACS (δ = 5 × 10 −6 ), BARON, HTS, and MHTS-TR approaches did not violate any constraints. In comparison with the original HTS, BARON, and CACS (δ = 5 × 10 −6 ) methods, the maximum profit obtained by the MHTS-TR approach was superior. Figure 7 plots the convergence graphs of the optimization results generated by the original HTS and MHTS-TR algorithms. It can be seen from Figure 7 that the proposed MHTS-TR approach showed a superior profit objective compared to the basic HTS method. In the convergence process, due to the appropriate trade-off between exploration and exploitation in our approach being performed by integrating different search mechanisms, the MHTS-TR algorithm converged to a feasible optimum very rapidly, meaning that the overall performance of the MHTS-TR approach was enhanced through the proposed modifications. In summary, the experimental results obtained by the MHTS-TR algorithm on this problem were better than those of the original HTS algorithm and the other competitors. Therefore, we can conclude that the MHTS-TR algorithm is applicable for solving real-world COPs.

Conclusions
Many real-world COPs are defined by complex mathematical equations with different constraints, and simply finding a feasible solution for such problems is not a straightforward task. Thus, to handle COPs efficiently, a novel approach with two search phases called MHTS-TR was proposed in this paper. The feasible search phase (the leader phase) ensured an intensified optimum in a relevant feasible region using the heat transfer search (HTS) algorithm, whereas the infeasible search phase (the follower phase) was used to introduce more diversification into the feasible search phase using the moving mechanism of the tandem running (TR) strategy.
To demonstrate the ability of the proposed MHTS-TR approach on handling different COPs, it was applied to a set of 24 constrained benchmark functions of CEC 2006, which involved different types of functions, such as, non-linear, linear, quadratic, polynomial, and cubic. Moreover, the simulation results were compared with those of the original HTS and other comparative algorithms. The obtained results showed that the proposed MHTS-TR approach markedly outperformed the other comparative algorithms. Finally, to validate the applicability of the MHTS-TR approach in solving real-world COPs, it was applied to handle a chemical COP, which was the simplified alkylation process. Moreover, the simulation results existing in the literature were introduced to be compared with the results obtained by our approach. The MHTS-TR algorithm demonstrated a superior result on the profit value compared to the original HTS method and some other competitors. Moreover, it showed competitive results in satisfying the constraints; therefore, we can conclude that the MHTS-TR approach is applicable for handling real-world COPs. The priority of each MHA is providing good solutions for specific problems, and therefore may not be as effective in providing solutions for some other problems. Thus, the considered application is intended for use on small-scale problems, and may not indicate towards the usefulness of the proposed approach for large-scale problems.
The future paths of research contain the application of the HTS-based optimization method for handling complex chemical engineering problems, such as multi-objective optimization problems involving a larger number of objectives and constraints. Moreover, to further check the robustness and suitability of the proposed method, it can be extended for solving other engineering processes, including single and multi-objective problems.