A new hybrid Lévy Quantum-behavior Butterfly Optimization Algorithm and its application in NL5 Muskingum model

1.   Introduction

Intelligent optimization algorithms are random search algorithms based on biological habits and natural phenomena. In recent years, intelligent optimization algorithms have been applied in many fields, such as medicine, engineering, and so on. Metanephritic algorithms can be classified into three main types: biological evolution, natural phenomena, and species' living habits. Biological evolution methods, like Genetic Algorithms (GA) ^[1] and Differential Evolution (DE) ^[2], are inspired by biological genetics, mutation, and evolution strategies. Natural phenomenon algorithms are based on the physical laws of nature, such as Sine Cosine Algorithm (SCA) ^[3] and Biogeography-Based Optimizer (BBO) ^[4]. Population life habits algorithms are inspired by the relationship between population individuals, including Particle Swarm Optimization (PSO) ^[5], Cuckoo Search Algorithm (CSA) ^[6], and Butterfly Optimization Algorithm (BOA) ^[7]. Scholars continue to develop new or improved algorithms despite the existence of many algorithms already because the $No \; Free \; Lunch \; Theory$ ^[8] argues logically that there is currently no optimization technique capable of solving all optimization problems. Additionally, when dealing with high-dimensional complex problems, metaheuristic algorithms tend to have slower convergence speeds and lower accuracy, making them susceptible to getting stuck in local optima ^[9]. Therefore, scientists persistently work on enhancing algorithms with superior performance to address experimental requirements such as modeling simulation ^[10], path planning ^[11], and robot frameworks ^[12].

BOA ^[7] is an intelligent optimization algorithm that imitates the foraging and courtship behavior of butterflies, and it has gained popularity among scientists and has been successfully applied in various fields. For example, Aygül et al. ^[13] applied BOA to calculate maximum power point tracking in PV systems under partial shadowing conditions. Arora and Anand ^[14] also used BOA to create learning automata. Li et al. ^[15] developed an updated version of BOA based on the cross-entropy technique to balance exploration and exploitation. However, its optimization performance significantly declines when faced with complex problems. This can be attributed to two main factors. First, there needs to be a balanced ratio between algorithmic exploration and search behavior exploitation throughout the search process. Second, the optimization process heavily relies on individual butterfly communication. When the optimal butterfly gets trapped in a local optimum, it attracts other butterflies, leading to premature convergence. In order to overcome these problems, scholars have taken some measures to improve the algorithm. Some scholars have rewritten the fragrance formula in the algorithm ^[16,17], which have some extent improved the convergence speed of the algorithm. However, the convergence accuracy and robustness still need to be improved. Assiri ^[18] attempted to incorporate the logistic mapping into the algorithm, but did not converge to the theoretical optimal value in function testing. Li et al. ^[19] applied quantum bits, quantum spiral gates, and quantum non-gates to the algorithm in a non-stage manner, improving its performance, but still having certain limitations when dealing high dimension problems. BOA also has been hybridized with other metaheuristics, such as Hybrid Butterfly Flower Pollination Algortihtm (HBFPA) ^[20] and Butterfly Artificial Bee Colony (BABC) ^[21].

Besides, the emergence of quantum computation theory has greatly influenced the field of computational science. Research focusing on incorporating quantum behavior into intelligent algorithms, such as the quantum genetic algorithm ^[22], the quantum PSO ^[23], and other new metaheuristic algorithms ^[24], has garnered significant interest. These advancements have successfully enhanced the efficiency and performance of these algorithms and helped scholars solve practical problems ^[25]. One notable application is the integration of quantum computation theory into optimization methods, with QPSO ^[26] being a popular example. The algorithm's ability to generate particles in a probabilistic manner throughout the solution space helps it avoid getting stuck in local optima, making researchers favor this, and often using QPSO as a tool for hybrid algorithms. Kumar et al. ^[27] proposed Cuckoo Search Adaptive Gaussian Quantum-behavior Particle Optimization (CSAGQPSO), which combines CSA and Adaptive Gaussian Quantum-behavior Particle Optimization (AGQPSO), a hybrid of two improved versions of QPSO, and applied it to ordinary differential equations. In CSAGQPSO, CSA and AGQPSO respectively undertake development and exploration tasks, which makes the algorithm simpler while having good performance. Mai et al. ^[28] combined QPSO and CSA, designing hybrid mechanisms to create a new hybrid Cuckoo-Quantum-behavior Particle Swarm Optimization (C-QPSO), successfully estimating the parameters of the Muskingum model. C-QPSO performs well and stands out in comparison with several improved algorithms, but the running cost of the algorithm is relatively high because there are many improvements in it. Wang et al. ^[29] introduced a novel QPSO variant by guiding it with the length of the potential well (LPW) and assisting it with Adam (abbreviated as Adam-LGQPSO), a combination of QPSO and Adam that enhances the algorithm's efficiency. The algorithm introduces the Adam algorithm in the initialization stage to improve the quality of the population and improve the optimization quality, but the effect is not very good in benchmark function experiments.

Despite its effectiveness, there is still room for improvement in terms of convergence accuracy. One tool that has shown promise in enhancing algorithm optimization performance is the concept of Lévy flight. Lévy flight is a random walk process that combines both local and global search strategies. In nature, many animals employ a mix of short and long distance searches to find food in unpredictable environments. Inspired by this, researchers have been attracted to the Lévy flight strategy and incorporated the idea of Lévy flight into swarm intelligence algorithms. Ling et al. ^[30] proposed an enhancement to the Wolf Optimization Algorithm method by integrating Lévy flight, which improved the global optimization capability of Lévy Wolf Optimization Algorithm. Zhong et al. ^[31] introduced the opposition-based learning Equilibrium Optimizer algorithm, which combines Lévy flight and evolutionary population dynamics to tackle high-dimensional global optimization problems. They have found that by allowing the algorithms to explore candidate solutions that are far away from the current optimal solution, the search space expands and the global search capability is enhanced, resulting in improved performance.

The hybrid algorithm resulting from the combination is more likely to incorporate the strengths of the two original algorithms while mitigating their weaknesses. Apart from HBFPA ^[20], BABC ^[21], CSAGQPSO ^[27], C-QPSO ^[28], and Adam-LGQPSO ^[29] mentioned previously, there are other successful instances as well ^[32,33]. By examining these successful cases, it is argued that the design of hybrid algorithms should emphasize the cooperative work between algorithms, with a key focus on inter-algorithm communication. The update mechanism of QPSO and BOA is simple and effective, which has a certain advantage in the design of hybrid algorithms. Second, in the update mechanisms of QPSO and BOA, the individual's position update is influenced by the $gbest$ point, promoting the cooperative work between the two algorithms and enabling them to leverage each other's strengths while minimizing their limitations. In addition, to maintain the computational speed of the algorithm, we divide the population into two equally large subgroups during the initialization phase, and design a hybrid mechanism to fuse the population during the update phase to further strengthen communication between algorithms.

Therefore, in this article, we propose a new hybrid algorithm that combines BOA and QPSO and uses Lévy flight, and we apply it to evaluate the parameters of the Muskingum model. BOA is an algorithm with a fast convergence rate, but it is prone to local optima, while QPSO is an algorithm with good global search ability, and the two can form functional complementarity. Specifically, the algorithm divides the population into two equally sized subgroups during the initialization phase, optimizes them through two algorithms, and then mixes the subgroups through a weight mixing mechanism during the optimization phase. In brief, the contributions of this paper are as follows:

● We propose a hybrid algorithm to estimate the parameters of the nonlinear five parameters Muskingum model, by combining both advantages of QPSO and BOA.

● We conduct extensive simulation experiments to evaluate our proposed algorithm, and the results show that our algorithm has better performance than some classic and recent offloading algorithms.

● We formulate the Muskinian model as parameter estimation problems of numerical solution to a differential equation, use new algorithms to evaluate model parameters, and improve the fit between evaluation values and observed values.

It should be pointed out that there has tremendous interest in developing hybrid intelligent algorithms to estimate the parameters of Muskingum models ^[34]. For example, Ouyang et al. ^[35] proposed a hybrid approach that combined PSO with the Nelder-Mead simplex (NMS) algorithm to optimize the parameters of the Muskingum model. Unlike traditional methods, the hybrid PSO-NMS approach does not require initial values for each parameter. Okkan and Kirdemir ^[36] combined the PSO algorithm with the Levenberg-Marquardt (LM) algorithm to estimate the parameters of a nonlinear Muskingum model with three parameters. Akbari and Hessami-Kermani ^[37] utilized the PSO-GA algorithm to optimize Muskingum parameters. It is important to note that the parameter optimization of the Muskingum model with four parameters (NL4) ^[28] is still in its early stages, and research on the Muskingum model with five parameters (NL5) ^[38] is even more limited. Thus, we will apply the proposed L-QBOA to estimate the parameters of the NL5 Muskingum model and hope to provide an effective new method for the study of the Muskingum model.

This article will describe L-QBOA in detail in Section 2, verify the performance of the algorithm through benchmark function testing in Section 3 and apply it to the NL5 Muskingum model, and finally summarize the full text in Section 4.

2.   A new hybrid algorithm

The main purpose of this paper is to present a new algorithm called L-QBOA, which is based on QPSO and BOA. The QPSO algorithm was modified and combined with BOA to create L-QBOA. To better understand L-QBOA, we will first provide a detailed description of BOA and QPSO.

2.1.   Butterfly Optimization Algorithm

Butterfly Optimization Algorithm (BOA) ^[7] is an intelligent algorithm that simulates the feeding behavior of butterflies. In this algorithm, there are $N$ butterflies.

where each butterfly $X_i$ represents a solution, $i = 1, \cdots, N$. The fitness of the butterfly is related to the solution vector of the butterfly. Butterflies can emit specific fragrances, communicate, and attract each other through the fragrances scattered in the air. Let $c$, $I$, and $a$ be the perception coefficient, the stimulus intensity, and the power exponent, respectively. The general value of the power exponent $a$ is 0.01. The equation of fragrance $f$ is given by

and then $c$ is given by

where $b$ is a constant with a value of 0.025, and $MaxIT$ is the maximum number of iterations. Each butterfly will fly randomly or move towards the one that emits a stronger fragrance. During the search process, butterflies have two modes of movement.

First, one of the butterflies emits the strongest fragrance and attracts them by being perceived by other butterflies. This is the global search stage for butterflies:

Second, when the scent in the air cannot be accurately perceived, they will be randomly attracted by a fragrance and move randomly. At this point, butterflies perform a local search:

Finally, butterflies will choose one of the modes for optimization in a probabilistic manner:

In above equations, $r$ is a random number $[0, 1]$, $p$ is usually set to 0.8, $gbest$ is the optimal solution in the current iteration, $f_i$ is the fragrance of the $i^{th}$ butterfly, and $X^t_j$ and $X^t_k$ respectively represent the solution vector of the $j^{th}$ and $k^{th}$ butterfly in the $t$ iterations.

2.2.   Modified Quantum-behavior Particle Swarm Optimization (MQPSO)

To improve the performance of QPSO ^[26], we add Lévy flight into QPSO by using the idea of reference ^[39]. Next, we give some details of the modification as follows.

In this algorithm, there are $N$ particles:

Each particle in the algorithm has corresponding local attractors $P_i$:

where $pbes{t_i}\left(t \right)$ is the particle's history optimal position, and $\varphi \left(t \right)$ is a random number in $[0, 1]$.

Here, Lévy flight, intensive short-range exploration and occasional long-range search, are used to expand the algorithm's search space and enrich the optimization methods of algorithms. Let $m(t)$ represent the mean of all particle's history optimal position and $u$ is a rand number in [0, 1]. When $u(t) \; >$ 0.5, particles will update their positions in Lévy flight:

where $\alpha$, called the contraction-expansion coefficient, decreases linearly from 1 to 0.5, and $L(\lambda)$ is Lévy flight:

where $r_a$ and $r_b$ are two standard normally distributed random numbers, $\beta$ is a constant equal to 1.5, and $\sigma$ is given by

$\Gamma \left(x\right)$ denotes the Gamma function. When $u(t) \; \leq$ 0.5, the $i^{th}$ particle's position is renewed by

2.3.   Lévy Quantum-behavior Butterfly Optimization Algorithm

The butterfly optimization algorithm exhibits a strong evolutionary mechanism. However, it also shares some limitations with other general metaheuristic algorithms. These include a slow rate of convergence, low precision, and a tendency to get trapped in local optima. To address these shortcomings, many hybrid algorithms have been developed to enhance performance by integrating the strengths of different algorithms. In this study, we propose a new algorithm, the Lévy Quantum-behavior Butterfly Optimization Algorithm (L-QBOA), by combining the above two algorithms and designing hybrid mechanisms. The Lévy quantum-behavior butterfly optimization algorithm will be explained in detail below.

In this method, the populations are updated by using BOA and MQPSO. First, the swarm is randomly divided into two equally sized sub-swarms $X$ and $Y$ during the initialization phase, and will not be re-divided thereafter. Then, swarm $X$ is updated by BOA and swarm $Y$ is updated using MQPSO, and their membership relationships remain unchanged in the subsequent iteration process. Later, L-QBOA performs the following mechanisms (weight hybrid mechanism, elite strategy and diversification mechanism) to continue the optimization work in the optimization phase.

2.3.1.   Weight hybrid mechanism

The algorithm divides the population into two subgroups during the initialization phase and optimizes them through two algorithms. To enhance collaboration between algorithms, a combination method is planned to be constructed during the optimization phase. Several combination ideas have been considered: first, sequentially performing different meta-heuristic algorithms overall; second, sequentially performing different meta-heuristic algorithms in each iteration ^[40]; and third, combining basic ideas of different meta-heuristic algorithms to design a new population updating strategy ^[41]. The second approach is deemed more suitable for BOA and QPSO, we argue.

Different from literature ^[40], weights are calculated for the two algorithms before mixing, rather than random mixing. Because the algorithm proposed in this article uses the basic BOA algorithm and QPSO with Lévy flight, there may be some differences in the quality of particles optimized through the two algorithms.

In each iteration, a new particle $Z$ is produced

where

or

and $f_X$ and $f_Y$ are the fitness of $X_i^t$ and $Y_i^t$, respectively. The expressions in Eq (2.14) or (2.15) can assign weights to $X$ in BOA and $Y$ in MQPSO to ensure that the new particle has a certain degree of randomness while appearing in a relatively good position. For example, use Eq (2.14) in minimizing optimization problems; if $Y^t_i$ is better, $f_Y$ will be smaller and $\eta _Y$ will be larger, allowing the information of particles with relatively good positions to be more fully retained in $Z$. If it is a maximization optimization problem, use Eq (2.15); if $Y^t_i$ is better, $f_Y$ will be larger and $\eta _Y$ will be larger.

Then, $Z$ is compared with $X$ and $Y$, keeping the better ones, and $pbest$ and $gbest$ in the MQPSO are updated. Here, the quality of $Z$ is not necessarily better than $X$ and $Y$, so in order to ensure the convergence of the algorithm, we can not blindly retain $Z$, so it is necessary to choose the best of $Z$, $X$ and $Y$ particles. In addition, the existence of $Z$ maintains the diversity of algorithm population to a certain extent, so this mechanism is multi-functional. This can lead the algorithm to develop in a positive direction.

2.3.2.   Elite strategy

Whether it is BOA or MQPSO, $gbest$ has always been utilized as a guiding factor in the algorithm's operation, thus the quality of its position greatly affects the algorithm's performance. To enhance the quality of $gbest$, we propose an elite strategy for its optimization. This strategy specifically targets the improvement of this position and generates a new random position near it, following a normal distribution

where $D$ represents the dimension of the problem, $range = \left({randg{e_1}, randg{e_2}, \ldots, randg{e_D}} \right)$, and $randge_d \; (d \in \left[ {1, D} \right])$ are the maximum distance values of the search space for each dimension, and $\xi$ is a larger positive number that can be adjusted according to different problems. Finally, we select the superior particle out of the two.

2.3.3.   Diversification mechanism

Population diversity is an important indicator that significantly affects algorithm performance. To maintain diversity, a diversification mechanism that mixes the mechanisms of literature ^[39] and literature ^[42] is employed in L-QBOA in the second half of the algorithm optimization work. The mechanism uses $NS$, a counter for the number of individual stagnation generations, recording the update of the best solution $gbest$. If $gbest$ remains unchanged in a given iteration, $NS$ is incremented by one; otherwise, it is reset to zero. When the value of $NS$ exceeds the set threshold $NS_{MAX}$, the algorithm determines that the population is overcrowded, resulting in a decrease in population diversity and falling into local optima. At this point, it is necessary to re-initialize ($\omega$%) of the poorer particles to maintain population diversity. What needs to be considered is that re-initializing too many particles will reduce the convergence speed of the algorithm, but re-initializing too few particles has little effect on maintaining the diversity of the algorithm. So, we need to set the two parameters according to the actual situation.

In the proposed L-QBOA described in this article, two independent populations are optimized by BOA and MQPSO separately. However, both algorithms jointly optimize one $gbest$, which represents the optimal solution from the two populations. This $gbest$ serves as a guiding factor for algorithm optimization, allowing the two independent algorithms to converge in the same direction. Real-time updates to the $gbest$ during the optimization process ensures algorithm convergence. Additionally, the $Z$ particles generated in the weight hybrid mechanism consist of $X$ and $Y$ particles. This inclusion helps maintain population diversity to a certain extent while retaining some information from $X$ and $Y$ particles. The algorithm establishes communication channels between the two algorithms using the $gbest$ and $Z$ particles, enabling each algorithm to function independently while maintaining their synergy.

The pseudo code of L-QBOA is shown in Algorithm 1. Meanwhile, the flow chart of the proposed L-QBOA algorithm is shown in Figure 1.

2.4.   L-QBOA performance analysis

In this part, this paper analyzes the computational complexity and memory efficiency of L-QBOA. These measure denote the time complexity and space complexity of the algorithm, respectively.

2.4.1.   Computation complexity

Computational complexity is a crucial metric in computer science that measures the difficulty of algorithms by assessing their efficiency and the computational resources they require. In general, the computational complexity of an algorithm tends to increase as the problem size grows, even if the algorithm remains the same. The complexity of algorithms plays a significant role in selecting appropriate algorithms or problem solutions. Hence, conducting a computational complexity analysis is essential to determine the feasibility and practicality of an algorithm for a given problem.

L-QBOA's computational complexity consists of four main components: population initialization, fitness value calculation, swarm location update, and the weight hybrid mechanism. We denote $N$ as the population size, $D$ as the dimension, and $T$ as the maximum number of iterations. During the swarm initialization phase, a random initialization generates a matrix, which requires a time complexity of $O(ND)$. Then, the fitness value of the population is calculated with a time complexity of $O(ND)$. When starting the iteration, the time complexity is proportional to $T$. In this phase, both the temporal complexity of fitness value calculation and location update are $O(ND)$. Finally, the weight hybrid mechanism requires a computational complexity of $O(ND)$. So, L-QBOA's computation procedure is described as

The computational complexity of BOA and QPSO can be easily calculated as $O(TND)$. So, the design of the new algorithm L-QBOA does not increase computational complexity.

2.4.2.   Memory efficiency

Spatial complexity is a metric that measures the memory resources used by algorithms. It quantifies the additional memory space needed during algorithm execution and how it increases with input size. When analyzing spatial complexity, we mainly consider the memory occupied by data structures, variables, pointers, and recursive calls. This analysis is essential in algorithm design and analysis as it helps evaluate resource consumption, choose suitable data structures, and optimize memory usage.

It is worth noting that the space complexity of the parameter setting of the algorithm is constant, that is, it is not affected by the complexity of the problem. The space complexity of L-QBOA is determined to be $O(N)$, where $N$ is the size of the problem. This complexity is mainly caused by the update of the population during iteration. The time complexity of the population is $O(N)$, so the space complexity of BOA and QPSO is $O(N)$.

In summary, based on the analysis, it can be concluded that L-QBOA and the two basic algorithms (BOA and QPSO) have similar space complexity characteristics, that is, linear space complexity. This means that the amount of additional memory space required by these algorithms grows proportionally to the problem size.

3.   Evaluation and experimental results

In this section, two class of numerical experiments are given to demonstrate the performance of our proposed hybrid algorithms. The first is the benchmark function experiment, and the other part is the NL5 Muskingum Model experiment.

3.1.   Experimental results of benchmark functions

All these algorithms are implemented by MATLAB R2016a and all runs are performed on a laptop with and an Intel(R) Core(TM) i5-7300HQ CPU @ 2.50 GHz, with 2*8 GB of RAM and with Windows 10 as the operating system. In the following experiments, the population size is set as 100, and the maximum number of iterations is 1000, and each operation is independently run 30 times. The parameter settings of BOA ^[7], MQPSO, QPSO + BOA, QSSA (2022) ^[43], HBA (2022) ^[44] and L-QBOA are shown in Table 1.

Before explaining the experiment, let us first explain an algorithm: QPSO + BOA. Compared to L-QBOA, this algorithm does not have the weight hybrid mechanism, elite strategy, and diversification mechanism. In QPSO + BOA, the $gbest$ in Eqs (2.8) and (2.4) is optimized by MQPSO and BOA jointly and BOA and MQPSO use $gbest$ as the information exchange channel. The purpose of sorting out such an algorithm is to more intuitively reflect the usefulness of our designed hybrid mechanism.

As is shown in Tables 2 and 3, we select 12 benchmark functions as the test set. F1–F6 are unimodal functions, and F7–F12 are multimodal functions. The characteristics of the unimodal functions and the multimodal functions are shown in these two tables, including the function name, expression of functions, dimension, the search area, and the theoretical optimal value. In order to illustrate the advantages of our proposed L-QBOA, we use different algorithms (BOA ^[7], MQPSO, QPSO + BOA, QSSA (2022) ^[43], HBA (2022) ^[44] and L-QBOA) to solve the minimum value of these benchmark functions. Here, in the following numerical experiments, we will analyze the performance of these algorithms from the following five aspects: convergency accuracy comparison, convergency speed comparison, robustness comparison, high dimension, and average CPU time.

3.1.1.   Convergency accuracy comparison

All statistical results viz, mean values (AVE), standard deviation (STD), maximum (MAX), and minimum (MIN) values of all benchmark functions and average CPU times are shown in Tables 4 and 5. The bold values in these tables represent the best values of the index in the function texts. Let us first examine the experimental data comparison between QPSO + BOA and BOA and MQPSO. The convergence accuracy of QPSO + BOA in F1–F8 is higher than that of BOA and MQPSO, which demonstrates that the combination of BOA and MQPSO can enhance the algorithm's performance. Even though MQPSO itself may be superior to BOA, their combination can further optimize each other's performance. When comparing L-QBOA and QPSO + BOA, it becomes evident that the performance of L-QBOA is significantly better. This proves that the hybrid mechanism can effectively enhance the algorithm's performance by balancing exploration and exploitation. Then we compare the performance of algorithms L-QBOA, QSSA, and HBA. Among the three algorithms, QSSA had the lowest performance, while HBA emerged as the winner in the F6 experiment. Interestingly, L-QBOA performed the best in all experiments except for F6. In F12, MQPSSO, QPSO + BOA, QSSA, and L-QBOA behaved the same. So, it can be concluded that L-QBOA has high convergence accuracy.

3.1.2.   Convergence speed comparison

The simulation graphs of algorithm iteration curves for 12 experiments are displayed in Figures 2 and 3. The convergence speeds of QPSO + BOA, MQPSO, and BOA were compared in 12 experiments. QPSO + BOA exhibited the fastest convergence speed, attributed to the cooperative effect of QPSO and BOA. It was observed that L-QBOA outperformed QPSO + BOA due to the enhanced convergence speed facilitated by the hybrid mechanism. QPSO + BOA has a certain degree of performance improvement compared to MQPSO and BOA, but due to the insufficient fine-grained optimization ability of algorithms, the algorithms cannot converge. When facing multimodal functions, BOA is prone to falling into local optima, making it difficult to converge to the optimal solution. In the initial eight experiments, QSSA exhibited slower convergence compared to HBA and L-QBOA. HBA demonstrated the fastest convergence in experiments F6 and F8, whereas L-QBOA outperformed in the remaining six experiments. The last four experiments showed comparable performance between the two algorithms. L-QBOA converges the fastest in most experiments. Consequently, it can be concluded that L-QBOA possesses a significant advantage in terms of convergence speed. (Notes: The algorithms appears to do not converge in most of the test functions, but in reality, the data output from MATLAB shows convergence (see Tables 4 and 5), which is a normal phenomenon ^{[30,31,45,46]}.)

3.1.3.   Robustness comparison

In order to compare the robustness of L-QBOA with BOA, MQPSO, QPSO + BOA and QSSA, we chose a search space of 30 dimensions with a maximum generation of 1000 for all test functions. The accuracy of all function experiments is 0, except for the Ackley function. Due to the particularity of the Ackley function, we set its corresponding experimental accuracy to 8.8818e-16. Each instance is calculated 30 times, and the relevant success rates and average iteration times of the five algorithms are shown in Table 6. The optimal values for each experiment are bold in this table. From the table, it can be seen that for unimodal function testing, in F1–F5, only L-QBOA can converge to the optimum with a success rate of 100%, while other algorithms have a success rate of 0, and only HBA can converge to the best in F6, but with a success rate of only 33.33%. In the case of multi-modal functions, QPSO + BOA has a higher success rate compared to BOA and MQPSO. Even when the success rate is the same, QPSO + BOA has a smaller average number of iterations compared to MQPSO and BOA. This is attributed to the synergistic effect of the two algorithms. The success rate of L-QBOA is the same as QPSO + BOA, but L-BOA has fewer iterations. It seems that the hybrid mechanism has played a role. Furthermore, for the multimodal functions, L-QBOA outperforms QSSA in terms of success rate and average number of iterations. For the success rate, L-QBOA has a 100% success rate of six, while HBA has only five. Although HBA requires an average of 1 less number of times in F9 experiments than L-QBOA, L-QBOA outperforms HBA in the remaining 5 experiments. Therefore, in terms of robustness, L-QBOA is superior.

3.1.4.   High dimensional testing

In order to test the performance of L-QBOA in solving high-dimensional optimization problems, the dimensions of the F1, F2, F6, F7, F11, and F12 were set to 100,500, and 1000, respectively for experiments, while other parameters settings remain unchanged, listed in Tables 7-9.

Table 4 illustrates that the performance of the F6 using the L-QBOA algorithm is inferior to that of HBA. However, when the dimension of function is increased to 100,500, and 1000, L-QBOA outperforms HBA. In experiments involving the F1 and F2, only the L-QBOA consistently computes the theoretical optimal value. As the dimension increases from 100 to 1000, all algorithms show a certain level of performance decline, except for L-QBOA which exhibits minimal decline.

When analyzing the data of the F7, F11, and F12 functions, combining Figures 4 and 5 can facilitate the comparison of algorithms. Figures 4 and 5 display the simulation convergence curves of three functions with dimensions of 500 and 1000, repectively. It is evident that L-QBOA algorithm demonstrates a clear advantage in computing F7. For F11, the performance of L-QBOA is on par with QSSA. However, in the case of F12, L-QBOA still has a slight advantage. Overall, L-QBOA has demonstrated good performance in computing high-dimensional optimization problems.

3.1.5.   Average CPU time comparison

For an intelligent algorithm, its intelligence is partly reflected in high optimization performance and fast running time. Part of our work is checking the running time of each algorithm. The average CPU time of all methods is shown in Table 10. It could be seen that QPSO + BOA runs in less time than MQPSO, indicating that dividing the population into two subgroups during the initialization phase can reduce algorithm running costs. While the hybrid mechanism may slightly increase runtime, in most experiments, the algorithm's execution time is still lower than that of the QSSA algorithm. Therefore, the design of L-QBOA focuses on enhancing algorithm performance while keeping operating costs manageable.

3.2.   Experimental results of the NL5 Muskingum model

To further verify the effectiveness and advantages of our presented algorithm, L-QBOA is applied to the parameter estimation problem of the NL5 Muskingum model

where $S(t)$, $I(t)$, and $O(t)$ represent channel storage, inflow, and outflow at time $t$, respectively; $K$ is the storage-time constant; $x$ is a weight factor; $\alpha _1$, $\alpha _2$ are flow parameters; and $m$ is an exponent parameter.

In order to estimate the parameters $K$, $x$, $\alpha _1$, $\alpha _2$, and $m$ in Eq (3.2), we use the following steps to construct an optimization problem:

Step 1: Assume $I(t_i)$ and $O(t_i)$ respectively represent the observed values of inflow and outflow, where $t_i = i/NT$ is the time nodal in a given time interval $[0, T]$ for $i = 0, 1, \ldots, N$. Otherwise, let $\hat O(t_i)$ be the calculation values of outflow, and the initial value is considered to be the same as the initial inflow: $\hat O(t_0) = O(t_0)$.

Step 2: For the time point $t_{i-1}$, we calculate $S(t)$ by using Eq (3.2).

Step 3: Use the backward difference formula to approximate the first-order derivative of Eq (3.1), and combining with Eq (3.2), we obtain

where $\Delta S({t_i}) = S({t_i}) - S({t_{i - 1}})$ and $\Delta {t_i} = {t_i} - {t_{i - 1}}$.

Step 4: Calculate the outflow in the next time

Step 5: Use the presented L-QBOA algorithm to solve the above optimization problem with SSQ as the objective function. The expression of SSQ will be given below (see Eq (3.5)).

In this paper, the performance of L-QBOA and the other parameter estimation procedures are evaluated by using the following measures:

1) SSQ: the sum of the square of the deviation.

2) SAD: sum of absolute values.

3) MARE: the mean absolute relative error.

4) PBIAS: the percent bias.

These four indicators can all be used to represent the fit degree of calculated data and observed data, and SSQ is generally used as the main reference indicator, which is also used as the objective function in this paper. The smaller the value, the higher the degree of fit between the evaluation value and the observed value.

This paper applied L-QBOA to three cases to verify the performance of our proposed method. The data calculated by SFLA-NMS and GA-GRG comes from reference ^[37]. L-QBOA also will be compared with C-QPSO, a new method once be applied to Muskingum model by Mai et al. ^[28] in 2023.

3.2.1.   Case 1

This is an example of parameter estimation for a nonlinear Muskingum model that Wilson ^[47] first considered. Table 11 shows the comparison of the SSQ, SAD, MARE and PBIAS for Wilson's data and optimal parameters obtained from L-QBOA and other methods (SFLA-NMS, GA-GRG and C-QPSO). The bold values in the table represent the optimal values of a metric calculated by algorithms. As demonstrated in Table 11, the SSQ = 5.442, SAD = 6.847, MARE = 0.253, and PBIAS = 0.29% calculated by L-QBOA in this case are all bolded in the table, which is the best value. In this case, L-QBOA estimates the parameters as $K = 0.837$, $x = 0.076$, $\alpha_1 = 0.695$, $\alpha_2 = 0.425$, and $m = 3.820$. As is noted, the computed optimal outflows by L-QBOA are better than those calculated by the other algorithms. The comparison between calculated values from L-QBOA and observed values is shown in Figure 6. It can be seen that L-QBOA performs significantly better than the other three algorithms.

3.2.2.   Case 2

The second example is a real example. This case study is a flood event that occurred in the Wye River, England. This flood event is a good example to validate the different processes of parameter estimation in the calibration step of the nonlinear Muskingum model. Table 12 displays the comparison of the best statistic values and corresponding parameters obtained from various algorithms. In this case, the PBIAS = 4.80% calculated by L-QBOA is the same as that calculated by C-QPSO, and the optimal value of PBIAS is 3.15% as calculated by FLA-NMS. However, the optimal values for SSQ, SAD and MARE are still provided by L-QBOA with 27804.91,602.723, and 2.584, respectively. {The corresponding parameters are $K$ = 0.466, $x$ = 0.687, $\alpha_1$ = 1.342, $\alpha_2$ = 1.049, and $m$ = 1.493.} The comparison between the outflow curve area and the observed values is shown in Figure 7.

3.2.3.   Case 3

The last example is the case of a flood with a multi-peak discharge hydrograph introduced by Viessman and Lewis in 2003 ^[49]. This particular case has a time interval $\Delta$ of 12 hours. The outflow hydrography exhibits two peaks occurring on the $10^{th}$ and $17^{th}$ days. The inflow ranges from 143 $cms$ to 1775.5 $cms$. The minimum and maximum outflow values are 118.4 and 1509.3 $cms$, respectively. Table 13 presents the SSD, SAD, MARE, and PBIAS values for the NL5 Muskingum model by using L-QBOA and other algorithms, respectively. The values in bold are the optimal values of the metric. In this case, C-QPSO provides an optimal value of PBIAS of 0.07%, and the calculated SSQ and SAD are better than that of SSFLA-NMS and GA-GRG. However, the optimal values for SSQ, SAD, and MARE are again calculated by L-QBOA as 65699.17,926.20, and 1.423, respectively. The corresponding coefficients are $K = 0.7954$, $x = 8.21\times10^{-8}$, $\alpha_1 = 3.364$, $\alpha_2 = 1.412$, and $m = 1.021$. The graph of outflow with the observed values is shown in Figure 8.

Tables 11–13 use SSQ as the main indicator, which is an objective function. In fact, other indicators derived from data may behave differently. When considering the impact of all four measures simultaneously, SFLA-NMS and C-QPSO face challenges. SFLA-NMS achieves a smaller PBIAS in Case 2, while C-QPSO achieves a smaller PBIAS in Case 3. However, both methods show a significant difference in SSQ from the optimal solution compared to L-QBOA. In conclusion, the parameter estimation performance of L-QBOA, as proposed in this paper for the NL5 Muskingum model, is as good or better than other methods.

4.   Conclusions

This article proposes L-QBOA, a hybrid algorithm that combines QPSO and BOA while introducing the Lévy flight mechanism. The algorithm divides the population into two subgroups during initialization and employs a hybrid mechanism during optimization to enhance collaboration. Experimental results on benchmark functions demonstrate that L-QBOA exhibits strong global optimization capabilities and performs well in high-dimensional problems. At the same time, it has been proven that using $gbest$ points to establish algorithm cooperative effects and a hybrid mechanism can effectively improve algorithm performance, which can serve as a reference for designing hybrid algorithms. To validate its practical applicability, L-QBOA was tested on the NL5 Muskingum model, yielding the best SSQ in all three cases. In conclusion, the proposed L-QBOA method shows advantages in the NL5 Muskingum model. In the future, L-QBOA can be considered for application to updated Muskingum models or other parameter estimation problems, numerical solutions of differential equations, and their inverse problems, path planning, or similar issues.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (12361087) and the Natural Science Foundation of Guangxi Province (2022GXNSFAA035618).

Conflict of interest

The authors declare there is no conflict of interest.