Improved Harris Hawks Optimization algorithm based on quantum correction and Nelder-Mead simplex method

1.   Introduction

Meta-heuristic (MH) optimization algorithms are extremely popular in engineering applications, because of superiority about simple concepts, easy implementation, no requirement for gradient information as well as capability of bypassing local optima, which can be treat as competitive solutions for various of problems covering different research fields, such as pattern recognition, target tracking, artificial intelligence and system control. By imitating social behaviors, biological principles or physical phenomena in nature, meta-heuristic algorithms can be categorized into 3 categories: evolution-based, physics-based and swarm-based approaches. A large number of scholars and researchers have successively proposed a variety of algorithms with specific strategies, for example, inspired by foraging behavior of Escherichia coli in human intestine, Passino ^[1] summarized the Bacterial Foraging Optimization (BFO) algorithm. Kennedy and Eberhart ^[2] proposed Particle Swarm Optimization (PSO) algorithm to find the optimal solution through cooperation and information sharing among individuals in the group. Differential Evolution (DE) algorithm based on differential simple mutation operation and one-to-one competitive survival strategy was presented by Storn and Price ^[3] in 1995. Moreover, Whale Optimization Algorithm (WOA) mimicking hunting behavior along with helix-shaped movement of humpbacks and Ant Colony Optimization (ACO) algorithm inspired by social intelligence of ants in finding the closest path from the nest to a source of food both attract considerable attention ^[4,5]. Besides, there are so many applications of meta-heuristic optimization algorithms for practical purposes. Yıldız wt al. ^[6] incorporated the chaotic maps in the elementary Lévy flight distribution and dubbed it as Chaotic Lévy Flight Distribution (CLFD) algorithm to address physical world engineering optimization problems. Wansasueb integrated Grey Wolf Optimizer (GWO), Genetic Algorithm (GA), Population Base Increment Learning (PBIL) and Water Cycle Algorithm (WCA) to form a new meta-heuristic named as Ensemble of Genetic algorithm (E-GGWP-W) to optimize design issue of composite wing ^{[7,8,9,10,11]}. Winyangkul et al. ^[12] presented multi-objective topology and sizing optimization of a morphing wing structure. Kumar et al. ^[13] enhanced 5 MH algorithms based on random migration search phase and simulated annealing-based selection for solving size and topology optimization of the trusses. Sparrow Search Algorithm (SSA) is employed to manage the operation of microgrid for minimizing either the total operating cost or the total emission ^[14]. Long et al. ^[15] proposed a novel hybrid seagull optimization algorithm (HSOA) based on cosine function and differential mutation strategy for estimating the parameters of photovoltaic models.

Harris Hawks Optimization (HHO) algorithm is a new kind of swarm intelligence algorithm proposed by Ali Asghar Heidari in 2019, which imitates the cooperative behavior of Harris hawks during predation ^[16]. The entire optimization process consists of three phases: exploration, exploitation and conversion. It enjoys the characteristics of simplicity in structure, concision in control parameters and excellent capability in global search. However, it suffers the shortcomings of low optimization accuracy, slow convergence speed and easy to fall into local optima. In order to settle the aforementioned issues, experts in different fields carried out a range of improvement method and put them into application. Tang et al. ^[17] promoted the convergence speed and accuracy of algorithm by introducing Tent Chaotic Map and elite hierarchy strategy into Chaos Elite Harris Hawks Optimization (CEHHO). Yin et al. ^[18] made use of nonlinear control parameter strategy and random inverse learning method to improve HHO (NOL-HHO). Attiya et al. ^[19] proposed an improved HHO algorithm (HHOSA) based on simulated annealing to solve problem of task scheduling in cloud computing. Ismael et al. ^[20] proposed Harris Hawks Optimization algorithm based on Opposition-Based Learning (HHOA-OBL) for hyperparameter estimation and feature selection. Qu et al. ^[21] used information exchange mechanism and nonlinear escape energy factors to improve population diversity and performance of algorithm. Ma et al. ^[22] mended fitness function through maximum likelihood estimation to get the solution for nonlinear equations about indoor positioning relate to difference of arrival time. Turabieh et al. ^[23] optimized HHO by controlling distribution of population and realized prediction of students' potential. ElSayed and Elattar ^[24] combined HHO with sequential quadratic programming to achieve optimal coordination for directional overcurrent relays of power. The Adaptive Harris Hawks Optimization (ADHHO) is proposed by Song for parameter identification of photovoltaic systems ^[25]. The persistent-trigonometric-differences mechanism and improved energy factor are helpful for balancing exploration and exploitation. With the help of comprehensive learning, equilibrium optimizer and terminal replacement mechanism, an improved algorithm named as comprehensive Learning Harris Hawks-Equilibrium Optimization (CLHHEO) is presented by Zhong for solving constrained optimization problems ^[26]. Hu et al. ^[27] used Specular Reflection Learning to improve HHO (HHOSRL) algorithm for assessing COVID-19. Dynamic Multi-Swarm Differential Learning Harris Hawks Optimizer (DMSDL-HHO) that divides population into many sub-swarms and introduces differential mutation operator candidate pool strategy is applied to the dispatch problem of hydropower stations ^[28]. Luo et al. ^[29] enhanced automatic epilepsy diagnosis method by using time-frequency analysis and Improved Harris Hawks Optimization (IHHO) with a hierarchical mechanism. Bardhan et al. ^[30] proposed ELM-IHHO by integrating the standard HHO algorithm with mutation-based search mechanism and Extreme Learning Machine for predicting soil compression index. Choi et al. ^[31] proposed an unsupervised intelligent system named as HHO-SVM for predicting the performance of a truck-haulage system using a combination of HHO and Support Vector Machine. Golafshani et al. ^[32] extended multi-layer neural network and radial basis function neural network by HHO to develop predictive models for the compressive strength of concretes containing supplementary cementitious materials.

As mentioned above, many scholars have applied variants of original HHO inspired by previous work in practical application, demonstrating the excellent performance and optimization capability of HHO. However, the original HHO still has potential for further improvement of convergence speed and solution accuracy. Moreover, there are serval defects belong to all or part of aforementioned algorithms. Firstly, the improved means for original HHO are not comprehensive, only focus on 2-3 aspects. There are opportunities to promote HHO in future. Secondly, significance of population's initial distribution is neglected, which directly affect the quality of subsequent optimization. Thirdly, there are some subjective factors in settings of parameters and constants, the reasons are not clearly stated, universality and robustness of the system may be reduced. Fourthly, according to a certain trigonometric function, individual takes movement in fixed pattern, it leads to the decrease of ergodicity and accuracy of search. Fifthly, the energy factor E used by the former studies does not accord with the law of physical consumption for animals. Finally, the influence between individuals in population is not taken into the consideration.

This paper takes efforts to improve original HHO and avoid above-mentioned problems as far as possible. First, change the generation mode of initial population. Chaos is a quite unique phenomenon caused by nonlinear effect, which has the characteristics of sensitivity to initial value, non-periodicity, long-term unpredictability and universality. So, chaos can exhibit the phenomenon of randomness and connection to the "random process" in system. Generating initial population by Hénon Chaotic Map based on conception of chaos can promote population's randomness and ergodicity in distribution. Second, improve the population diversity and local search capability of HHO. Quantum computing uses quantum bits code to store data, for the same size of population, compared with other methods, quantum code can contain much more information. So, introducing quantum correction in exploitation phase under certain condition can make use of features of quantum computing to promote diversity and increase the opportunities of taking precise movement of individual, under the premise of size of population is no changed. Third, enhance global search ability of HHO. Modified Nelder-Mead simplex method is widely used to solve parameter estimation and similar statistical problems, which requires only one or two function evaluations per iteration and quickly produces satisfactory results in first few iterations. So, modified Nelder-Mead simplex method is used to boost search performance and breadth in exploration phase, empower individual to move in wide range as needed in the search space, reducing the probability of falling into local optima. At the same time, it will not excessively increase computational complexity. Fourth, group communication factor belong to BFO can express relationship between individuals effectively. Attraction and repellent power are introduced into the search process as parameters can optimize function evaluation, and then, population regeneration will be more accurate. Final, when animal is running or jumping, with energy consumption, muscle cells will produce lactic acid, then muscles show fatigue. At this time, the speed or jumping distance will be declined, animal need to stop and have a rest. Integrate the conception of cycle about "run—consume physical strength—stop to take rest—restore energy" into update of escape energy factor E can better balance exploitation and exploration of HHO in accordant with biological laws.

The improved HHO studied in this paper is named as Quantum Correction Harris Hawks Optimization (QC-HHO) algorithm. It is tested by 10 classic benchmark functions and 30 CEC2014 benchmark functions, experimental results are recorded. In order to provide accurate and reliable conclusion, Wilcoxon rank-sum test with 5% degree is carefully performed. Compared with other optimizers, it can be generalized that the QC-HHO outperforms others in most cases. Moreover, for the application of gas leakage source localization, QC-HHO can provide excellent experimental result, search accuracy and speed are satisfactory.

2.   Introduction of Harris Hawks Optimization algorithm

Harris Hawks Optimization algorithm is a meta-heuristic algorithm derived from Harris hawk's predation action. This section briefly describes the basic principles and mathematical models of it.

2.1.   Exploration phase

The exploration phase of HHO algorithm is global search process. When hawks find the target prey in the air, all individuals coordinate their action to find a favorable position around the prey and form siege.

In initial state, the individual Hawk appears in a certain position belong to search space according to the principle of randomness, then gradually moves to the optimal solution. Set q as a random number located between 0 and 1, when q < 0.5, each individual will take movement by reference the position of other ones and prey. As q ≥ 0.5, the Harris hawk will stay in a tree that is in the scope of activity of population. The mathematical model can be described as follows:

In Eq (2.1), X(t) is position of population after the t^th round of iteration. X_rand is location of a random individual in population. X_Prey is position of the optimal solution, X_m is the average position of population. r₁, r₂, r₃ and r₄ are the random numbers in the interval of (0, 1). Ub and lb are the upper and lower bounds of search space respectively, N is population size.

2.2.   Exploitation phase

The exploitation phase is a process of local search. When the prey is surrounded, the hawks will take attack. The HHO algorithm simulates escape behavior of prey and hunting strategy of hawks through different combinations of escape energy factor E and random numbers r ∈ (0, 1). r is escape opportunity for prey. r < 0.5 means that the escape is successful, but r ≥ 0.5 denotes an unsuccessful escape.

1) Soft besiege

As |E| ≥ 0.5 and r ≥ 0.5, prey is energetic and try to escape by jump, but eventually it is caught, a mathematical formula is as follows:

In Eq (2.3), ΔX(t) represents the distance between optimal solution and current individual after the t^th time of iteration, r₅ is a random number between 0 and 1. J is jumping distance of prey during running for life.

2) Hard besiege

As |E| < 0.5 and r ≥ 0.5, prey is lack of physical strength and captured directly, formula is as follows:

3) Soft besiege with progressive rapid dives

As |E| ≥ 0.5 and r < 0.5, prey is energetic and has chance to escape, hawk will make a more intelligent soft encirclement, the implementation is as follows:

In Eq (2.8), D is the dimension of problem, S is a D-dimensional random vector, and LF is the Levy flight function as follows:

In Eqs (2.10) and (2.11), ru and rv are random numbers between 0 and 1, β is a constant numbered as 1.5.

4) Hard besiege with progressive rapid dives

As |E| < 0.5 and r < 0.5, prey has insufficient physical energy, but it still has a chance to escape. In order to reduce the average distance from prey, hawks form a new hard encirclement, the strategy is as follows:

2.3.   Conversion phase between exploration and exploitation

The HHO algorithm controls the conversion between global search and local search through the escape energy factor E, which is defined as follows:

In Eq (2.15), T is the maximum rounds of iterations, E₀ as (-1, 1) is the initial value of the energy in iterations. When |E| ≥ 1, the HHO algorithm enter into the exploration phase, which represents individuals in population move towards the prey within whole search space. When |E| < 1, search will convert into the exploitation phase.

3.   Description of QC-HHO algorithm

In order to overcome the shortcomings of original HHO, as well as improving the searching performance, this section will propose the optimization methods in 5 aspects include Hénon Chaotic Map, quantum correction, modified Nelder-Mead simplex method, group communication factor and amended escape energy factor E.

3.1.   Hénon chaotic map

This section will give the definition and mathematical model of the Hénon Chaotic Map, as well as compare the distribution of initial population generated by random number and Hénon Chaotic Map.

High-quality initial population can help to improve the accuracy and speed of convergence. If the distribution of initial population is far from the target solution, it is impossible to effectively solve the problem ^[33]. Therefore, in the process of initialization, the initial values should be distributed as evenly as possible in the solution space in order to get better results. The original HHO algorithm is initialized by random numbers, which cannot guarantee the diversity and ergodicity of population.

As complex behavior of nonlinear systems, chaos exhibits phenomenon that is similar to randomness. This feature can be used to improve the performance of the algorithm ^[34]. Hénon Chaotic Map is proposed by Michel Hénon in 1976, who is a French mathematician ^[35]. It is discrete-time dynamic system which can generate chaotic phenomenon in 2-Dimensional space, the mathematical expression is as follows:

In Eq (3.1), when a takes 1.4, b takes 0.3, system will enter into chaotic state. The initial state of Harris hawk population is as follows:

In Eq (3.2), N is scale of population, D is dimension of problem. Figure 1 shows the initial population generated by two methods. (a) is the population generated by random numbers, and (b) is generated by Hénon Chaotic Map. It can be seen that (b) is better than (a) in both ergodicity and randomness.

3.2.   Quantum correction

This section uses quantum bits code and quantum rotation gate to improve the population diversity and local search capability of HHO. In 1981, quantum computing was proposed by Richard Feynman who is an American theoretical physicist. The basic unit of information in storage is called as a quantum bit. Each quantum bit may represent 1 or 0, even the superimposed state of 0 and 1, which can be expressed by a linear superposition of two orthogonal basis vectors:

α and β are probability amplitudes, which represent the linear probability of individual states $\left|\text{0}⟩\right.$ and $\left|\text{1}⟩\right.$ respectively and satisfy following conditions:

One quantum bit can be expressed as $\left[\begin{array}{c}\text{α}\\ \text{β}\end{array}\right]$. n quantum bits can describe 2ⁿ state at the same time. Therefore, for the same size of population, compared with other methods, quantum code can contain more information. Because evaluating optimization algorithm through fitness function have to use a certain value, it is necessary to perform a measurement for each quantum bit that will collapse into a certain value. The mathematical formula is as follows:

For the purpose of enhancing the population diversity and search accuracy, the position X of individual in population is converted into a quantum bit that can be updated by quantum rotation gate, the process is as follows:

1) Obtain parameter P that is the ratio of individual position X to search space.

For example, position of individual x_i ∈ X is (7, 9), ub and lb are 15 and 6 respectively. So, the value of P is as follow:

2) Convert a decimal integer ${I}{{P}}_{{i}}^{{j}}$ = ⌊${{P}}_{{i}}^{{j}}{}{×}{}{1000}$⌋, i = {1, 2, ..., N}, j = {1, 2, ..., D} into a 10-bit binary number ${{BP}}_{{i}}^{{j}}{}{ = }{}\left(\begin{array}{cccc}{{b}}_{{1}} & {{b}}_{{2}} & {⋯} & {{b}}_{{10}}\end{array}\right)$, N is population size and D is dimension of problem. So Eq (3.7) can deduce as follow:

3) Take a random number r ∈ (0, 1), calculating α_k and β_k according to Eqs (3.4) and (3.5) respectively, ${{BP}}_{{i}}^{{j}}$ will be converted into a quantum bit sequence that contain 10 quantum bits as follows:

The entire population ${X} = \{{{x}}_{{i}}^{{j}}{}\}$ can be expressed as $\left\{{{QB}}_{{i}}^{{j}}\right\}$, i ∈ {1, 2, ..., N}, j ∈ {1, 2, ..., D}. Based on assumption of r = 0.4, Eq (3.9) can be expressed as follow:

4) Quantum rotation gate is used to change the distribution domain of α and β, and its form is as follows:

θ in Eq (3.13) is rotation angle, and the update process of quantum bit is as follows:

For example, set θ as 0.2, after one time of rotation, ${{QB}}_{{i}}^{{j}}$ in Eqs (3.11) and (3.12) will be as follow:

5) According to the quantum state updated by Eqs (3.14) and (3.5), the quantum bit ${{QB}}_{{i}}^{{j}}$ is collapsed into individual position ${X}{}{'}{}{ = }{}{{}\{{{x}}_{{i}}^{{j}}{}}\}{'}$, the fitness of function ${f}{(}{X}{'}{)}$ is calculated, as well as update X_Prey(t). So, ${{QB}}_{{i}}^{{j}}$ in Eqs (3.15) and (3.16) can be convert to BP_i as follow:

6) If the count of rotations does not reach the maximum m, return to Step 4).

Since the measurement, collapse and rotation of the quantum bit will increase the time consumption, in order to balance the accuracy and efficiency of optimization, the quantum rotation operation is only performed in the situations that adopting soft besiege strategies or hard besiege strategies, meanwhile, the curve convergence state is as smooth state. Compared with HHO, this method gives the search agent more opportunities to take delicate movement, so it is more likely to get the better results.

3.3.   Modified Nelder-Mead simplex method

Nelder-Mead technique is proposed by John Nelder and Mead Roger in 1965, it is a kind of heuristic search method with simplex concept that can converge to a non-stable point for the problem that can be solved by an alternative method. It is usually used for the nonlinear optimization problem of which derivative is unknown and issues for finding the minimum or maximum value of objective function in a multi-dimensional space ^[36].

The NM method constructs an initial simplex containing a given point, and then the point with the worst value of function is replaced by way of reflection, expansion and contraction. If three methods mentioned above are all failed, contraction reduce the radius of the simplex to small enough. The modified NM method is introduced in global search phase to enhance search ability of HHO. The steps are as follows:

1) Set X_NM-Optimal is the optimal point in population, X_{NM-Suboptimal} is the suboptimal point, and X_NM-Worst and X_NM-Subworst are the worst point and the point that is next to the worst point respectively.

2) Generate a group of new positions X_NM-Elite ^[17] with parameters of X_NM-Optimal, X_{NM-Suboptimal,} X_NM-Subworst and X_NM-Worst, the formula is as follows:

3) Let X_NM = {X_NM-Elite^T, X_NM-Optimal, X_{NM-Suboptimal}, X_NM-Worst, X_NM-Subworst} denote the list of points in the simplex.

4) Sort the points X_NM in the simplex from lowest function value to highest one. Update X_NM-Optimal, X_NM-Subworst, X_NM-Worst at each step in the iteration.

5) Set reflected point X_NMr and expanded point X_NMe as follows:

In Eqs (3.19) and (3.20), X_m is the average position of population in Eq (2.2).

6) If f(X_NMr) < f(X_NM-Optimal), X_NM-Worst is replaced according to follow rules:

a. If f(X_NMe) < f(X_NMr), accept X_NMe, named as "Expansion".

b. Otherwise, accept X_NMr, named as "Reflection <2> ".

7) If f(X_NM-Optimal) ≤ f(X_NMr) < f(X_NM-Subworst), X_NM-Worst is replaced by X_NMr, named as "Reflection <1> ",

8) If f(X_NM-Subworst) ≤ f(X_NMr), perform "Contraction" by follow rules:

a. If f(X_NMr) < f(X_NM-Worst), calculate the contracted point X_NMc1 as follows:

If f(X_NMc1) < f(X_NMr), X_NM-Worst is replaced by X_NMc1, named as "Outside contraction", otherwise, continue with Step 9).

b. If f(X_NMr) ≥ f(X_NM-Worst), set the contracted point X_NMc2 as follows:

If f(X_NMc2) < f(X_NM-Worst), accept X_NMc2, named as "Inside contraction", otherwise, continue with Step 9). The individual X_worst corresponding to X_NM-Worst in population is updated.

9) Update list of X_NM:

Calculate f(v(i)), i = 2, 3, ..., 8. The simplex at the next iteration is {X_NM(1), v(2), ..., v(8)}, named as "Shrink".

10) Go back to Step 4), until the stop condition is met.

11) If f(X_NM-Optimal) < f(X_Prey(t)), set Fitness(t) = f(X_NM-Optimal), and update X_Prey(t) as X_NM-Optimal.

Take an example, set cost function as:

For convenience of understanding and expressing the procedure of NM method, the simplex in this example contains only 3 points, X_NM = {X_NM-Optimal, X_NM-Worst, X_NM-Subworst}, initial values of them are [-2.9723, 0.1679], [-1.9723, 0.1679] and [-2.9723, 1.1679]. From Step 1) to Step 3) are skipped, after 19 times of iteration, data are shown in Table 1.

As shown in Table 1, each iteration possesses 2 lines of data, the top line is the status of f(X_NM) after execution of Step 4), the bottom one is the result of Step 9).

1) At the 2^nd iteration, fit for Step 6a), so action of "Expansion" is executed, X_NM-Worst is changed and f(X_NM-Worst) is 3.54 that is shown in the bottom line of 2^nd iteration.

2) At the 3^rd iteration, fit for Step 7), action of "Reflection <1> " is executed, f(X_NM-Worst) is changed to 3.54.

3) At the 4^th iteration, fit for Step 9), X_NM is updated.

4) At the 5^th iteration, fit for Step 8b), action of "Inside contraction" is executed, f(X_NM-Worst) will be equal to f(X_NMc2).

5) At the 12^th iteration, fit for Step 8a), action of "Outside contraction" is executed, f(X_NM-Worst) is changed as the values of f(X_NMc1).

6) At the 19^th iteration, fit for Step 6b), action of "Reflection <2> " is executed, X_NM-Worst is equal to X_NMr.

3.4.   Group communication factor

Harris hawk is a kind of animals hunting in group. When they encircle prey, there are 2 or 3 hawks stop in one tree at most, it means that position of each hawk is affected by another one. Original HHO does not take the relationship between individuals into consideration. The attraction and repulsion between creatures are introduced into HHO. The attractive force makes hawks to move toward to same prey. The repulsive force keeps hawks at a certain distance and not gather in a single position. The group communication factor gcf is added to the fitness of cost function in order to generate a new fitness for comparison and optimization. The group communication factor gcf of the i^th search agent is defined as follows:

In Eq (3.25), N is population size, d_attract is the attraction released by hawk, and ω_attract indicates the diffusivity of attraction. h_repel and ω_repel are the influence and width of repellent respectively. The variable vd is the distance parameter between different variables of the i^th search agent, which is defined as follows:

In Eq (3.26), D is the dimension of search agent. ${\text{X}}_{\text{i}}^{\text{m}}$ is the position of the m^th variable of the i^th search agent.

In Eq (3.27), N is population size, Fitness_gcf is fitness value of function that integrated with group communication factor after iteration. Therefore, introduction of gcf improves comprehensiveness and completeness of fitness.

3.5.   Escape energy factor E based on bio-energy consumption

In HHO algorithm, escape energy factor E reflects the search ability for the optimal solution, determines the conversion between global search and local search as well as affects the strategy adopted in exploitation phase. As shown in Eq (2.15), the value of E is depended on the random number E₀, it appears as linear decreasing trend including many times of fluctuation during optimization process, which can hardly describe actual change tendency of creatures' energy. Some scholars improved factor E from various of aspects, such as the m-HHO algorithm proposed by Gupta et al. ^[37] set nonlinear energy factor E as follows in order to obtain more opportunities for exploitation.

In Eq (3.28), t is count of iteration, T is the maximum number of iterations, mi is the nonlinear modulation index, and E_start and E_end are the initial and final energy parameter values respectively.

Although m-HHO get better results in exploration phase, the universality and robustness of algorithm are reduced because of abandoning random parameter. So HHO and m-HHO are both not in line with the energy consumption and recovery of prey during escape or actual situation of predator and prey. Prey's physical fitness decreases rapidly with the increase of running distance in escape. It can recover to a certain extent through a short rest. QC-HHO separates the change process of E in HHO into several groups, each group contains corresponding data of E in iteration, that is as following formula:

In Eq (3.29), r is a random number in the interval of [50, 100] that means the size of each group. In Eq (3.30), j determines the frequency of physical recovery of prey in escape, namely, how many times of rest the prey takes for. E₁_Part_j is the escape energy factor of each group in whole search process. In Eq (3.32), E₀ is in the interval of [-1, 1] that is as same as it of HHO in Eq (2.15). In Eq (2.5) relate to HHO, jump distance of prey is irrelevant to energy. However, according to Eq (3.33), the more energy the prey has, the longer distance it can jump.

As shown in Figure 2, (a) is the change trend of escape energy factor E of original HHO, which is relevant linear, (b) is the curve of E in m-HHO algorithm, and (c) is E of QC-HHO algorithm. It can be seen that the fluctuation of E in each group is getting smaller and smaller. The energy of prey is recovery after short break when it is exhausted. The physical consumption of prey in the process of being chased is decreasing, the upper limit of physical will be lower and lower as time goes on. That is conform to the feature of bio-energy.

4.   Pseudo code and spatiotemporal complexity analysis of algorithm

4.1.   Pseudo code and flow chart of QC-HHO

This section describes the application of methods mentioned in former sections as form of pseudo code, the corresponding flow chart is shown in Figure 3.

Input: N is population size, T is maximum number of iterations

Output: Position of prey and corresponding fitness value

Execute Eqs (3.1) and (3.2) to initialize the population X_i(i = 1, 2, ..., N).

Initialize the parameters of Nelder-Mead simplex method

While (stopping condition is not met) do

  Check state for convergence curve

  Execute Eq (3.27) to calculate the fitness values integrating with the group communication factor gcf

  for (each hawk X_i) do

  Execute Eq (3.25) to calculate the group communication factor gcf

   Execute Eq (3.32) to update the escape energy factor E

  Execute Eq (3.33) to update the jump strength J

  if1 the state of convergence curve is rapid convergence and iteration is in later stage

    set E = 1 to force the search enter into exploration

  endif1

  if2 (|E| ≥ 1) then exploration phase

    if3 convergence curve is rapid convergence

      Execute Eq (2.1) to update the location

    else

      Execute Eqs (3.18)–(3.23) to use Nelder-Mead simplex method

    endif3

  endif2

  if4 (|E| < 1) then exploitation phase

    if5 (r ≥ 0.5 and |E| ≥ 0.5) then soft besiege

      Execute Eq (2.3) to update the location

    endif5

    if6 (r ≥ 0.5 and |E| < 0.5) then hard besiege

      Execute Eq (2.6) to update the location

    endif6

    if7 (r ≥ 0.5 and it is possible to drop into local optima)

      Execute Eqs (3.6), (3.10), (3.13) and (3.14) to execute quantum correction

    endif7

    if8 (r < 0.5 and |E| ≥ 0.5) then soft besiege with progressive rapid dives

      Execute Eq (2.9) to update the location

    endif8

    if9 (r < 0.5 and |E| < 0.5) then hard besiege with progressive rapid dives

      Execute Eq (2.14) to update the location

    endif9

  endif4

endfor

endwhile

Return X_Prey

4.2.   Analysis of time complexity

Set the population size of the QC-HHO is N and iteration rounds is T, D is dimension of problem, the time complexity of algorithm according to the execution steps of the algorithm is analyzed.

Step 1, initialization of hawk population, the time complexity is O(2N × D).

Step 2, initialize the parameters of Nelder-Mead simplex method, the time complexity is O(1).

Step 3, calculate the convergence trend, the time complexity is O(1).

Step 4, calculate the function fitness, and the time complexity is O(1).

Step 5, perform group behavior.

Step 5.1, calculate the group communication factor gcf, the time complexity is O(N).

Step 5.2, calculate the escape energy factor E and jump strength J, the time complexity is O(2).

Step 5.3, judge whether to force search enter into exploration according to the convergence trend, make comparison once, make calculation once, the time complexity is O(2).

Step 5.4, confirm to enter into exploration or exploitation, make comparison once, the time complexity is O(1).

Step 5.5, in exploration phase, compare twice, move once, update parameters of Nelder-Mead simplex method, and the time complexity is O(4).

Step 5.6, in the exploitation phase, compare once to determine the attack strategy, move once, the time complexity is O(1).

Step 5.7, in soft besiege, move once, execute quantum correction, the time complexity is O(20 × D + 12).

Step 5.8, in hard besiege, move once, execute quantum correction, the time complexity is O(20 × D + 12).

Step 5.9, in soft and hard besiege with progressive rapid dives, move once, the time complexity is O(1).

Step 5 will repeat for N times, so its time complexity is O((10 + N + 20 × D + 12) × N).

After the above steps, the time complexity of QC-HHO after T iterations is O(T ×((10 + N + 20 × D + 12) × N)).

4.3.   Analysis of spatial complexity

The spatial complexity is evaluating indicator of through storage space used by algorithm and analyzed according to the procedures. Set the scale of population is N, the number of iterations is T, and the dimension of function is D. X[N][D] stores the position of population, X_Prey[1][D] saves the optimized variable, Fitness keeps the optimal result for fitness of function, gcf is the relationship between individuals in population, NM_X[7][1] stores the parameters of Nelder-Mead simplex method, NM_X_Fitness[2][1] is fitness of function belong to NM simplex method, P[N][D] keeps ratio of individual position X to search space, IP[N][D] stores decimal integer of P, BP[10] are binary data of individual, QB[N][D][2][10] stores quantum bits of population, X'[N][D] stores the position of population after quantum correction. Therefore, the spatial complexity of QC-HHO algorithm is: O(N × D × 20) + 4 × O(N × D) + 2 × O(D) + O(21).

5.   Test and analysis

This section makes use of 10 classic benchmark functions and 30 CEC2014 benchmark functions to test the performance of QC-HHO on the platform of intel i7-4790, 8GB DRAM and MatLab 2012b, and presents the experimental results and analysis.

5.1.   Comparison on classic benchmark functions

This section is to evaluate the performance of QC-HHO, taking comparison between QC-HHO and BFO ^[1], PSO ^[2], DE ^[3], HHO ^[16] as well as other 4 improved HHO algorithms including CEHHO ^[17], ADHHO ^[25], CLHHEO ^[26], DMSDL-HHO ^[28] on classic benchmark functions.

Classic benchmark functions are grouped into three categories: F1−F4 are uni-modal functions that used to evaluate the capability of exploitation. F5−F8 are multi-modal functions, F9 and F10 are fixed-dimensional multi-modal functions. F5−F10 have more than one local optima, the number of variables increases exponentially, which are used to evaluate the ability of exploration and avoiding the local optima. The figures, names, variable number, variable boundary and theoretical optimal solution of functions are shown in Table 2.

In detail, for the tests on classic benchmark functions, the termination criterion is the maximum number of iterations T is 500, the number of populations is 30, the rounds of runs on each function is 30, the maximum number m of quantum rotations is 3, means and standard deviations are collected as support for evaluation and shown in Tables 3 and 4.

As shown in Table 3, the performance of QC-HHO is significantly better than the other 4 classic optimization algorithms, especially for F1, F2, F3, F6, the theoretical optimal solution can be found. For uni-modal functions F1−F4, compared with the other algorithms, the means and standard deviations of QC-HHO are not only good at final result, but are more stable. In terms of multi-modal functions and fixed-dimensional multi-modal functions, the mean of F5, F7, F9, QC-HHO is similar to HHO and superior to BFO and DE. For F10, QC-HHO obtains satisfactory results, and have a slight advantage over HHO.

As shown in Table 4, theoretical optimal value of F1, F2, F3 are found by QC-HHO and DMSDL -HHO that completely suppress other 3 improved HHO algorithms about best value and mean, they are both ranked first. stability of DMSDL-HHO in F3 is the best among all 5 algorithms. For F4, QC-HHO is ranked second and slightly behind CLHHEO, but it is more stable. For F5, the result of QC-HHO is not only outstanding, but also is more stable substantially. For F6 and F9, 5 optimizers get theoretical optimal values without exception. F7, results of 5 algorithms are identical, but they are not theoretical optimal value. That is means that model of original HHO is not suitable for this function. For F8 and F10, QC-HHO get ranking the second and third. QC-HHO, it can be seen that the capability of solving multi-modal problems is not as good as this of uni-modal ones.

The convergence curves of classic benchmark functions are shown in Figures 4−6. The blue solid line with rectangle is the QC-HHO, the purple solid line is CEHHO, the green solid line is ADHHO, the black dashed line is DMSDL-HHO and the red dot dashed line is CLHHEO.

Figure 4 is convergence curves of uni-modal benchmark functions. Among the 5 improved HHO algorithms, QC-HHO has the best final optimal result. Convergence speed of F1 is the fastest, better than the other 4 algorithms, and finds the theoretical optimal solution. For F2 and F3, although QC-HHO is slower in convergence speed in the early phase of search, speed is significantly faster in the later phase, and a better result is found through quantum rotation correction in the final phase. For F4, QC-HHO fell into the local optima at least 2 times, but all of them successfully jumped out and keep searching the search space to find good solutions.

Figure 5 is comparison of multi-modal benchmark functions. For F5 and F6, QC-HHO has a slightly slower convergence rate than ADHHO and CLHHEO, and achieves a jump out of the local optima in initial phase of F5. For F7, under the premise of identical final result, the convergence speed of QC-HHO is the fastest among 5 algorithms. F8, QC-HHO's convergence speed is not outstanding, and the final result is worse than this of ADHHO.

Figure 6 shows convergence curves of 5 improved HHO in fixed-dimensional multi-modal functions F9 and F10. The final optimization results are at the same order of magnitude, but QC-HHO jumps out of local optima for several times. For F9, the final results of all optimizes are as same as theoretical optimal value, the only differences between them are convergence rate. QC-HHO takes one time of jumping out of local optima to get better result.

5.2.   Comparison on CEC2014 benchmark functions

This section compares QC-HHO to CEHHO ^[17], ADHHO ^[25], CLHHEO ^[26] and DMSDL-HHO ^[28] on 30 CEC2014 benchmark functions for further tests of QC-HHO's performance.

For CEC2014 test, population size N is set to 100, the search range is in [−100,100]^D, for each function, each algorithm independently run 50 times to get experimental results that composed of mean and standard deviation values on 30, 50, and 100-dimensional CEC2014 benchmark functions in Tables 5, 7 and 9, respectively. Moreover, for convenience of making analysis, ranking statistics for each category are recorded in Tables 6, 8 and 10.

As shown in Tables 5 and 6, the results obtained in 30-dimensional CEC2014 benchmark functions show that QC-HHO obtains 15 times ranking the first, 5 times ranking the second, 4 times ranking the third, 5 times ranking the fourth and 1 time ranking the sixth, and no ranking the fifth, its average ranking is 2.1 which is less than this of any other optimizers. DMSDL-HHO and ADHHO compute the average rankings are 2.17 and 2.30 respectively, which are also relatively competitive. The reason why QC-HHO can be the top one is that it is extremely good at Unimodal Functions and Simple Multimodal Functions. QC-HHO is outstanding in average ranking of these two categories, it is means quantum correction mechanism and modified Nelder-Mead simplex method guarantee QC-HHO can find better solution in exploration and exploitation phase by precise or long-distance movement, even if other optimizers are incapable of get better results. Also, the modified escape energy factor E makes the conversion between global and local search at right time. Distribution of initial population generated by Hénon Chaotic Map makes search more accurate. Uniform and moderate distance between individuals is good for group communication factor play a more effective role in the process of population regeneration. On the other side, intermediate level optimization ability result in QC-HHO is surpassed by DMSDL-HHO and ADHHO in categories of Hybrid Functions and Composition Functions.

As shown in Tables 7 and 8, on the 50-dimensional CEC2014 benchmark functions, top one is also QC-HHO that obtains 14 times ranking the first, 5 times ranking the second, 5 times ranking the third, 4 times ranking the fourth and 2 time ranking the sixth, and no ranking the fifth, its average ranking is 3.23 that is more than this of 30-dimension to some extent. No.2 to No.4 positions are also DMSDL-HHO, ADHHO and CLHHEO. The competitiveness of HHO and CEHHO are relatively weak in comparison to the others.

As shown in Tables 9 and 10, on the 100-dimensional CEC2014 benchmark functions, QC-HHO obtains 15 times ranking the first, 6 times ranking the second, 3 times ranking the third, 5 times ranking the fourth and 1 time ranking the fifth, and no ranking the sixth, its average ranking is 2.03 that is better than this of 30 and 50-dimension. Increase of dimension makes QC-HHO take the advantage of distribution of initial population effectively. Results of Composition Functions are obviously promoted. Higher dimension makes taking long-distance movement of individuals by Nelder-Mead simplex method easier, be instrumental in bypass the local optima. As the result, the final solution is better than lower-dimensional population.

According to the analysis aforementioned, QC-HHO has good performance in solution accuracy and more reliable scalability than other competitors, be able to get the most times on ranking the first, and its average rankings are always the highest in category of all functions. Specifically, QC-HHO achieves the highest average ranking on categories of Unimodal Functions and Simple Multimodal Function, however, for Hybrid Functions and Composition Functions, QC-HHO is incapable of achieving the top one, just position in the second or the third, that proves that QC-HHO has a robust exploration and exploitation capability for finding an optimal solution, at the same time, it need to be enhanced in future study.

5.3.   Analysis of time consumption

This section takes analysis about efficiency and time consumption of QC-HHO on classic benchmark functions and CEC2014 benchmark functions. The results are shown in Tables 11 and 12 that correspond with the testing data in Tables 3−10.

In order to assess the computational complexity of QC-HHO, the mean values of the runtime over all runs of each algorithm on classic benchmark functions and 30-dimensional CEC2014 benchmark function are recorded in Tables 11 and 12. The numbers in the brackets behind the time values are the rank of computational complexity for every function.

The runtime taken by 6 kinds of optimizers to find the solutions for F1–F10 of classic benchmark functions are listed in Table 11. As per results in the table, QC-HHO shows a reasonably fast and competitive performance in finding the best solutions compared to most of optimizers for uni-modal cases. It takes 4 ranking the first, 2 ranking the second, 3 ranking the third, 1 ranking the fourth, 1 rank fifth, which is belong to moderate levels. F4 is an exception, all 5 improved HHO take more time than HHO, but QC-HHO has relatively good performance. For multi-modal and fixed-dimensional multi-modal cases F5−F10, the search efficiency of QC-HHO is slightly lower than this of HHO, excluding for F7. This result imputes to execution of quantum correction in exploitation phase. Quantum correction can take more exact and accurate search to bypass optima, but it is a kind of time-consuming operation. On the other side, QC-HHO is superior to most of other 4 improved HHO algorithms in most cases. As a whole, although QC-HHO takes less computational cost only for part of uni-modal functions, it performs faster than HHO, CEHHO, DMSDL-HHO, ADHHO, CLHHEO. Because of quantum correction mechanism and modified Nelder-Mead simplex method, time consumption of QC-HHO on multi-modal functions is very high.

Each algorithm's average runtime results on all functions, and the ranking results according to the average runtime are shown in the last two rows. The results in Table 12. Because procedures of CEHHO, DMSDL-HHO, ADHHO, CLHHEO and QC-HHO are all more complicated than original HHO, the runtime of most cases are longer than these of HHO. Average runtime of DMSDL-HHO and QC-HHO are 5.7069 and 5.8127 respectively, the rank of QC-HHO is 3 that is lower than this of DMSDL-HHO. However, it is just lower in 2%, at the same time, the mean for QC-HHO on the 30-dimensional CEC2014 benchmark functions is higher than this DMSDL-HHO in 8%. It is worth of exchange 2% of the time cost to 8% of the performance improvement. Furthermore, QC-HHO spend more time than other optimizers in Hybrid Functions and Composition Functions, QC-HHO just rank third and fourth. Because the functions in these two categories are much more complicated than Unimodal Functions and Simple Multimodal Functions, quantum correction as a kind of time-consuming operation will be executed more frequently. Also, in order to bypass the local optima, individual need to move in long distance by modified Nelder-Mead simplex method. These two optimization mechanisms are both increase the computational complexity.

5.4.   Result of Wilcoxon rank-sum test

In order to investigate the significant differences between the results of proposed QC-HHO versus other optimizers and provide more accurate and reliable conclusions, Wilcoxon rank-sum test with 5% degree is carefully performed. Tables 13 and 14 shows the attained p-values of the Wilcoxon rank-sum test with 5% significance.

5.5.   Application for gas leakage source localization

This section is applied to problem of gas leakage source localization. The results of QC-HHO are compared to various modified optimizers proposed in previous studies.

This problem can be described mathematically as follows:

Consider: z = [x, y, q]

Variable range:

Where

In Eq (5.1), ${\text{C}}_{\text{SensorData}}^{\text{i}}$ is the data of gas concentration collected by the i^th sensor. ${\text{C}}_{\text{Try}\text{Data}}^{\text{i}}$ is the data that is predicted by Gaussian puff model in Eqs (5.2) and (5.3). f(z) means the deviation of concentration value between predicted by Gaussian puff model and measured by sensor. The smaller f(z), the closer to actual gas leakage position. ξ is Eq (5.2) is random error simulated by white Gaussian noise with standard deviation of ${\text{σ}}_{\text{wge}}^{\text{2}}$. X_Sensor and Y_Sensor are the location of sensor.

The experimental site is set as Figure 7, there are 9 sensors set in it. Assuming the coordinate of gas leakage sources is (2, 3), gas release rate q is 30ml/min. the mean velocity of wind speed is 12.50 cm/s. The standard deviation σ_wge of measurement error is set to 0.5, the final positioning results are shown in Table 15.

Although the index of optimal cost of QC-HHO is not outstanding in all optimizer, the deviations of position and release rate are excellent and smaller than these of others. Good results are obtained at an acceptable cost of time.

6.   Conclusions

This paper briefly explains the principles of classic HHO and proposes QC-HHO that includes 5 aspects of improvement: optimizing the distribution of initial populations by Hénon Chaotic Map, promoting local search and population diversity by quantum correction, enhancing global search performance by Nelder-Mead simplex method, describing the relationship between individuals by group communication factor, improving the selection of strategies in exploitation phase as well as occasion of conversion between global and local search by modified escape energy factor E based on biological energy consumption. A test was conducted on 10 classic benchmark functions and 30 CEC2014 benchmark functions to analyze exploration, exploitation, capability of jumping out of local optima, and convergence feature of the proposed algorithm. The data and convergence curve show that QC-HHO is competitive with other optimization algorithms in finding the theoretical optimal solution, it improves the efficiency and robustness of original HHO. On the other hand, the performance of the multi-modal benchmark function is not as good as that of uni-modal benchmark function, which needs to be further studied.

For application of indoor gas leakage source localization through wireless sensor networks, the experimental data shows that the accuracy of position and gas release rate are excellent, but index of time consumption is just in the middle level. As the whole, the performance corresponds to the result of testing of benchmark function. Limited by objective conditions and project requirements, QC-HHO is not testing in outdoor environment. Also, for multiple leakage source, if distance between leakage sources is too close and the diffusion concentration fields are overlapped, there will be large deviation in the locating results. Furthermore, quantum correction is helpful for jumping local optima, but it is a kind of extremely time-consuming operation, so the time when to execute it should be more accurate.

Acknowledgments

This article is funded by: Research on Key Issues of Gas Leakage Source Location in Compressed Sensing Sensor Network in Three-dimensional Complex Flow Field Project of Tianjin Natural Science Foundation (20JCYBJC00320).

Conflict of interest

All authors declare there is no conflicts of interest.