An Enhanced Multioperator Runge–Kutta Algorithm for Optimizing Complex Water Engineering Problems

Abstract

Water engineering problems are typically nonlinear, multivariable, and multimodal optimization problems. Accurate water engineering problem optimization helps predict these systems’ performance. This paper proposes a novel optimization algorithm named enhanced multioperator Runge–Kutta optimization (EMRUN) to accurately solve different types of water engineering problems. The EMRUN’s novelty is focused mainly on enhancing the exploration stage, utilizing the Runge–Kutta search mechanism (RK-SM), the covariance matrix adaptation evolution strategy (CMA-ES) techniques, and improving the exploitation stage by using the enhanced solution quality (IESQ) and sequential quadratic programming (SQP) methods. In addition to that, adaptive parameters were included to improve the stability of these two stages. The superior performance of EMRUN is initially tested against a set of CEC-17 benchmark functions. Afterward, the proposed algorithm extracts parameters from an eight-parameter Muskingum model. Finally, the EMRUM is applied to a practical hydropower multireservoir system. The experimental findings show that EMRUN performs much better than advanced optimization approaches. Furthermore, the EMRUN has demonstrated the ability to converge up to 99.99% of the global solution. According to the findings, the suggested method is a competitive algorithm that should be considered in optimizing water engineering problems.

1. Introduction

Extensive research on optimization methods in various areas of water engineering systems, such as reservoir operations, flood management [1,2,3,4], groundwater management [5,6,7], quality management of water resources [8,9,10,11], water distribution systems [12,13,14], and sedimentation [15,16,17], has attracted the attention of researchers to improve optimization algorithms to solve complex water engineering systems. Metaheuristic optimization methods (MOMs) have piqued the interest of many academics in recent years and have been extensively used for various real-world purposes. In general, most MOMs are motivated by biological or physical behaviors. For example, the particle swarm optimization (PSO) [18], Harris hawks optimization (HHO) [19], and the grey wolf optimizer (GWO) [20] are based on biological behaviors, while heat transfer search (HTS) [21], differential evolution (DE) [22], gradient-based optimization (GBO) [23], the weighted mean of vectors (INFO) [24], and the gravity search algorithm (GSA) [25] are motivated by certain basic physical or mathematical concepts. MOMs provide several benefits over typical gradient-based algorithms, including better global search ability, easier implementation, and no need for gradient information [26,27,28,29,30].

The approaches described above have been used to handle a variety of water engineering optimization issues [31,32,33]. A hybrid harmony search (HHS) approach was presented by Karahan et al. for estimating nonlinear Muskingum model (NMM) parameters [34]. The Broyden, Fletcher, Goldfarb, and Shanno (BFGS) algorithm is used as a local search method to speed up the HHS process. From the findings, it can be concluded that the suggested approaches can be used to confidently predict the optimum parameter values of the NMM. Bozorg-Haddad et al. (2015) developed a hybrid shuffling frog leaping algorithm (SFLA) and the Nelder–Mead simplex (NMS) for optimizing a novel Muskingum model with a seven-parameter set of constraints [35]. A comparison of the results of the methodology of this study with other models found that their methodology performed significantly better. Modified honey bee mating optimization (MHBMO) developed by Niazkar and Afzali was used to determine the parameters used for identifying Muskingum models (MUMs) [1]. The MHBMO was efficiently used in two case studies to estimate the parameters of the MUM. It is shown by the obtained findings that the proposed hybrid technique can be used effectively for estimating the stringent nonlinear Muskingum (NLMs) parameters. To extract the optimal parameters from an NMM with four parameters, Farahani et al. used the Kidney algorithm (KA) [36]. According to their findings, the KA algorithm is an excellent choice for optimizing the NMM with high precision. The teaching-learning-based optimization (TLBO) method was used by Bozorg-Haddad et al. to estimate an NMM with three parameters [37]. The findings of the Nash–Sutcliffe efficiency (NSE) measure demonstrated that the suggested technique performed very well in calculating the Muskingum parameters with high precision.

Jothiprakash and Arunkumar used the chaotic technique with the DE and the genetic algorithm (GA) to optimize hydropower output from a reservoir system [38]. Chaotic DE (CDE) outperformed other strategies in terms of yearly power generation. Hosseini-Moghari et al. used the imperialist competitive algorithm (ICA) and cuckoo search algorithm (CS) to optimize the Karoun-4 reservoir operation in Iran [39]. The findings showed that CS could become closer to the global optimum. A novel hybrid algorithm has been devised by combining the GA with the krill algorithm, according to Ehteram et al., which is intended to optimize the operation of multireservoir systems [40]. The suggested hybrid method was tested by using a reservoir from Iran and two benchmark issues (i.e., 4-reservoir and 10-reservoir systems). According to the study results, the suggested hybrid algorithm outperformed the other approaches. To identify the optimal operating guidelines for multireservoir systems, Ahmadianfar et al. proposed a multistrategy DE and PSO (MS-DEPSO) method that combined DE and PSO [26]. According to the findings, the suggested algorithm outperforms the competition on its capacity to derive the higher number of optimal operating principles for reservoir systems. Using a multimechanism ensemble interior search algorithm (MEISA), Abdollahi and Ahmadianfar attempted to solve a variety of challenges in the optimal management of reservoir systems [41]. The findings revealed that the MEISA performs very well at establishing the best policies for the effective management of multireservoir systems.

The methods described above have effectively addressed numerous optimization issues (for example, water engineering problems). Although they have produced virtually acceptable results, they have some benefits and drawbacks. The most significant benefit of these methods is their flexibility, in that they can be modified to incorporate different nonlinear systems, including a specific set of constraints and fitness functions. Several significant drawbacks include difficulty dealing with the randomized pattern of algorithm parameters, slow convergence, and an inability to discriminate global solutions.

Because of the nature of the complexity of water engineering optimization issues, the development of new, efficient techniques suited for tackling these optimization problems is required. Most MOMs use different powerful operators and suitable adjustments in the control parameters to achieve a balance between global and local searches. Despite this, researchers have not yet discovered a way to effectively exchange information across MAs’ global and local skills. In 2021, Ahmadianfar et al. invented the Runge–Kutta (RK) optimization (RUN) [27], a population-based approach free of metaphors.

When compared with other MOMs, RUN has the following distinguishing characteristics: (i) a single factor aids RUN in making a seamless transition from the global to local phase, (ii) the RK method offers better global and local searches in solution spaces, and (iii) the enhanced solution quality (ESQ) boosts the quality of solutions. Although present MOMs can address specific optimization issues, they cannot tackle all optimization issues. When the exploration efficiency of the RUN improves, the exploitation capacity decreases, and vice versa. Therefore, obtaining a correct balance between these two phases through the optimization technique is a concern, especially in certain complex cases that need an enhanced optimization method. As a result, the current study intends to create an enhanced multioperator RUN (EMRUN) algorithm that enables the optimal search to match the exploration and exploitation phases.

Furthermore, an updated version of ESQ is used to eliminate placement in specific local solutions. Finally, the covariance matrix adaptation evolution strategy (CMA-ES) [42] and sequential quadratic programming (SQP) [43] are used as an effective and robust local search to assist the EMRUN technique in accelerating the convergence and accuracy of solutions. The EMRUN algorithm is evaluated using 29 test functions and real-world water engineering challenges to determine its effectiveness.

The following outlines the remainder of this paper’s organization: Section 2 describes the basic RUN and the proposed EMRUN algorithm. Section 3 evaluates the performance of the EMRUN on CEC-17 benchmark functions. The proposed algorithm is assessed by extracting optimal parameters for the Muskingum model in Section 4. Section 5 investigates the performance of EMRUN for optimizing a hydropower multireservoir problem, followed by Section 6, which closes this research article.

2. Enhanced Multioperator RUN Algorithm

2.1. The RUN Algorithm

The Runge–Kutta technique is the foundation of the RUN algorithm [27]. Two operators have been used in this algorithm: the Runge–Kutta search method (RK-SM) and the ESQ. The primary processing steps in the RUN algorithm are outlined in the subsequent subsections.

2.1.1. Updating the Location of Solutions

Based on the RK approach, the RUN method uses the RK-SM to change the position of the candidate solution (xn) at each iteration. The main formulation for the RK-SM is defined as follows:

$$x_{new}=\{\begin{array}{c}\left(x_{m1}+v.AC.r.x_m\right)+AC.\left(RK-SM\right)+\mu .randn.(x_{m2}-x_{m1}) if rand<0.5\\ \left(x_{m2}+v.AC.r.x_{m2}\right)+AC.\left(RK-SM\right)+\mu .randn.\left(x_{b1}-x_{b2}\right) otherwise \end{array}$$

(1)

in which $x_{m1}$ and $x_{m2}$ are given by

$$x_{m1}=\beta .x_{avg}+\left(1-\beta \right).x_{b1}$$

(1a)

$$x_{m2}=\beta .x_{best}+\left(1-\beta \right).x_{lbest}$$

(1b)

and

$$AC=2.(0.5-rand)\times \alpha$$

(1c)

$$\alpha =10\times exp\left(-12.rand.\left(\frac{Iter}{MaxIter}\right)\right)$$

(1d)

where $v$ denotes an integer number with a random amount of 1 or −1; $x_{b1}$ and $x_{b2}$ express two randomly selected solutions of [1, Np]; Np denotes the number of populations; $r$ denotes a random number from the interval [0, 2]; $AC$ is an adaptive coefficient, and $\mu$ indicates a random number. $x_{best}$ indicates the best-so-far solution, and $x_{lbest}$, is the best position calculated at each iteration. $RK-SM$ is expressed as follows:

$$RK-SM=\frac{1}{6}.\left(X_{RK}\right).\Delta x$$

(1e)

$$X_{RK}=w_1+2\times w_2+2\times w_3+w_4$$

(1f)

$$w_1=\frac{1}{2\Delta x}\left(rand.x_{\mathrm{wst}}-\eta .x_{\mathrm{bst}}\right)$$

(1g)

$$\eta =round\left(1+rand\right).\left(1-rand\right)$$

(1h)

$$w_2=\frac{1}{2\Delta x}\left(rand.\left(x_{\mathrm{wst}}+e_1.k_1.\Delta x\right)-(\eta .x_{\mathrm{bs}}+e_2.k_1.\Delta x\right))$$

(1i)

$$w_3=\frac{1}{2\Delta x}(rand.(x_{\mathrm{wst}}+e_1.\left(\frac{1}{2}{w}_2\right).\Delta x)-(\eta .x_{\mathrm{bst}}+e_2.\left(\frac{1}{2}{w}_2\right).\Delta x))$$

(1j)

$$w_4=\frac{1}{2\Delta x}\left(rand.\left(x_{\mathrm{wst}}+e_1.w_3.\Delta x\right)-(\eta .x_{\mathrm{bst}}+e_2.w_3.\Delta x\right))$$

(1k)

where $x_{\mathrm{bst}}$ and $x_{\mathrm{wst}}$ indicate the best and worst solutions, respectively, and $e_1$ and $e_2$ indicate two random numbers in the interval [0, 1]. $\Delta x$ is formulated as follows:

$$\Delta x=2.rand.\left|Q\right|$$

(1l)

$$Q=rand.\left[\left(x_{bst}-rand.x_{avg}\right)+\rho \right]$$

(1m)

$$\rho =rand.\left(x_n-rand.\left(U-L\right)\right).\exp (-4.\frac{Iter}{MaxIter})$$

(1n)

where $Q$ expresses a step size; $L$ and $U$ express the lower and upper limitations of search space for each problem; $x_{avg}$ is the average of solutions; $Iter$ denotes the iteration number; and $MaxIter$ denotes the maximum number of iterations. In this paper, $x_{bst}$ and $x_{wst}$ are determined as follows:

$$\begin{gathered} if f\left(x_n\right)<f\left(x_{bs,n}\right) \\ {x}_{bst}=x_n \\ {x}_{wst}=x_{bs,n} \\ else \\ {x}_{bst}=x_{bs,n} \\ {x}_{wst}=x_m \\ end \end{gathered}$$

(1o)

where $x_{bs,n}$ indicates the best solution achieved among three random solutions ($x_{b1}$, $x_{b2}$, and $x_{b3}$).

2.1.2. ESQ Operator

RUN uses an effective operator, ESQ, to increase solutions’ quality while avoiding local solutions. The RUN algorithm uses the ESQ operator and generates the solution $U$ using the following scheme:

$$\begin{gathered} if rand<0.5 \\ U=\{\begin{array}{c}{x}_2+\lambda .\rho .|\left(x_2-x_{avg}\right)+randn| if \lambda <1 \\ \left(x_2-x_{avg}\right)+r.\delta .|\left(2.rand.x_2-x_{avg}\right)+randn| otherwise\end{array} \\ end \end{gathered}$$

(2)

$$\delta =rand\left(0, 2\right).exp\left(-p.\left(\frac{Iter}{MaxIter}\right)\right)$$

(2a)

$$x_{avg}=\frac{x_{b1}+x_{b2}+x_{b3}}{3}$$

(2b)

$$x_2=\theta \times x_{avg}+\left(1-\theta \right)\times x_{best}$$

(2c)

where $\theta$ denotes a random number in the interval [0, 1], $p$ denotes a random number with the amount of 5$\times rand$, and $r$ expresses a random integer number, equal to 1, 0, or −1.

Because it is possible that solution $U$ does not have a superior objective function than solution $x_n$ (i.e., $f\left(U\right)>f(x_n)$), solution $x_3$ is produced to give the problem a second opportunity to develop a viable solution, which is expressed as follows:

$$\begin{gathered} ifrand< \lambda \\ {x}_3=\left(U-rand.U\right)+AC.(rand.X_{RK}+\left(2.rand.x_{best}-U\right)) \\ end \end{gathered}$$

(3)

2.2. The Proposed EMRUN Algorithm

In this paper, the enhanced EMRUN is developed to optimize the water engineering problems. The suggested method uses a rank-based mechanism to pick the most effective approach for updating the solutions, an adaptive mechanism to update the algorithm parameters, and an improved ESQ (IESQ) operator to increase the efficacy of solutions and accelerate the convergence process. Two powerful local search algorithms, CMA-ES and SQP, are eventually implemented. The following sections describe the basic structure of EMRUN.

2.2.1. Search Mechanism of EMRUN

The suggested EMRUN uses an efficient search mechanism (ESM) to investigate the potential locations inside the viable space. The ESM selects the exploration and exploitation phases at each iteration by applying an effective rank-based algorithm. The ESM can be represented as follows:

$$V_1=\{\begin{array}{c}\left(x_{m1}+v.A{C}^{\prime }.x_m\right)+A{C}^{\prime }.\left(RK-SM+0.5.\sigma .(x_{m2}-x_{m1}\right)) if rand>{\psi }_n\\ \left(x_{m2}+v.A{C}^{\prime }.x_{m2}\right)+A{C}^{\prime }.\left(RK-SM+0.5.\sigma .\left(x_{b1}-x_{b2}\right)\right) otherwise \end{array}$$

(4)

in which

$$\sigma =0.5+0.1.randn$$

(4a)

$$A{C}^{\prime }=4.(0.5-rand)\times {\alpha }^{\prime }$$

(4b)

$${\alpha }^{\prime }=2\times exp\left(-8.rand.\left(\frac{Iter}{MaxIter}\right)\right)$$

(4c)

and

$$\psi \left(I{d}_n\right)=1-\cos \left(P\left(I{d}_n\right).\frac{\pi }{2}\right)$$

(4d)

with

$$P\left(I{d}_n\right)={\left[\frac{\left[Np-n\right]}{Np}\right]}^2$$

(4e)

where $\psi$ indicates the rank of each solution, $P$ denotes the overall rank of each solution, and $I{d}_n$ indicates the solution’s index after sorting on the basis of the objective function.

Additionally, the crossover operator is used to increase the variety of the population. The solution produced by crossover ($x_1$) can be expressed as follows:

$$x_{1,k}=\{\begin{array}{c}{V}_{1,k} if rand<p{c}_{1,n} \mathrm{or} k=k_{rnd} \\ x_{n,k} otherwise \end{array}$$

(5)

in which

$$p{c}_{1,n}={\varphi }_{pc1,n}+0.1.randn$$

(5a)

where $p{c}_{1,n}$ expresses the crossover probability factor, $k_{rnd}$ indicates an integer number with randomness nature in the interval [1, D], and ${\varphi }_{pc1,n}$ indicates a number with a value equal to 0.5 at the initial iteration. ${\varphi }_{pc1,n}$ is calculated according to Equation (5b):

$${\varphi }_{pc1,n}=\omega .{\varphi }_{pc1,n}+\left(1-\omega \right).\frac{{{\displaystyle \sum }}_{i=1}^{size\left(Qp{c}_1\right)}Qp{c}_{1,i}}{size\left(Qp{c}_1\right)}$$

(5b)

where $\omega$ denotes a number with an amount equal to 0.9 and $Qp{c}_1$ indicates all successful amounts of $p{c}_{1,n}$ during each iteration. During each iteration of the EMRUN, the best amount of $p{c}_{1,n}$ for each population can be employed, which can then be utilized to find relevant places in the viable space for each individual.

2.2.2. Improved Enhanced Solution Quality (IESQ)

EMRUN leverages the IESQ operator to increase the solution quality and prevent local solutions. The suggested EMRUN uses two differential vectors that are based on two random solutions and the average solution. As a result, the EMRUN algorithm uses the IESQ to construct the solution $V_2$ by using the following procedure:

$$\begin{gathered} if rand<{\left(1-Iter/MaxIter\right)}^2 or It<0.1.MaxIter \\ {V}_2=\{\begin{array}{c}{x}_{b1}+A{C}^{\prime }.\left(x_{avg}-x_n\right). B if rand<\tau \\ x_n+g_n.\left(x_{BS}-x_n\right)+0.5.g_n.\left(x_{btr,r1}-x_{por,r2}\right) otherwise \end{array} \\ end \end{gathered}$$

(6a)

in which

$$x_{BS}=L.x_{btr,r1}+\left(1-L\right).x_{avg}$$

(6b)

$$B=rand\left(1,D\right)<\tau$$

$$\tau =0.5+0.1.randn$$

(7)

where $B$ is a logic factor to determine which component of the differential vector $\left(x_{avg}-x_n\right)$ can be used to create the solution $V_2$. $x_{BS}$ is a solution based on $x_{btr,r1 }$ or $x_{avg}$. $L$ denotes a random number with a value equal to 0 or 1. If $rand<0.5$, L is equal to 1; otherwise, it is equal to 0. $x_{por,r2}$ and $x_{btr,r1}$ denote the powerless and better random solutions (where $r1\ne r2\ne n$). To define these solutions, the population is first split into two groups with equal members (Np/2). The first group consists of solutions with a higher rank, based on their objective functions, while the second group consists of solutions with a lower rank. As a result, the solutions $x_{btr,r1}$ and $x_{por,r2}$ are chosen randomly from the first and second groups, respectively. To determine the $g_n$ as an adaptive coefficient, it needs to calculate the solution $x_2$. Therefore, the solution $x_2$ is first calculated below.

Furthermore, the crossover operator is used to produce the new solution that is based on Equation (8):

$$x_{2,k}=\{\begin{array}{c}{V}_{2,k} if rand<p{c}_{2,n} \mathrm{or} k=k_{rnd} \\ x_{n,k} otherwise \end{array}$$

(8)

in which

$$p{c}_{2,n}={\varphi }_{pc2,n}+0.1.randn$$

(8a)

where $p{c}_{2,n}$ indicates the crossover probability factor and ${\varphi }_{pc2,n}$ indicates a number with a value equal to 0.5 at the initial iteration. ${\varphi }_{pc2,n}$ is calculated according to Equation (8b):

$${\varphi }_{pc2,n}=\omega .{\varphi }_{pc2,n}+\left(1-\omega \right).\frac{{{\displaystyle \sum }}_{i=1}^{size\left(Qp{c}_2\right)}Qp{c}_{2,i}}{size\left(Qp{c}_2\right)}$$

(8b)

where $Qp{c}_2$ indicates all successful $p{c}_{2,n}$ during each iteration. The adaptive coefficient $g_n$ is expressed as

$$g_n={\varphi }_{g_n,n}+0.1.\tan \left(0.5.\pi .randn\right)$$

(8c)

with

$${\varphi }_{g_n,n}=\omega .{\varphi }_{g_n,n}+\left(1-\omega \right).\frac{{{\displaystyle \sum }}_{i=1}^{size\left(Q_{g_n}\right)}{\phi }^{\prime }{Q}_{g_n,i}{}^2}{{{\displaystyle \sum }}_{i=1}^{size\left(Q_{g_n}\right)}{\phi }^{\prime }{Q}_{g_n,i}}$$

(8d)

$${\phi }^{\prime }=\frac{\left|f\left(x_2\right)-f\left(x_n\right)\right|}{{{\displaystyle \sum }}_{i=1}^{size\left(Q_{g_n}\right)}\left|f\left(x_{2,i}\right)-f\left(x_{n,i}\right)\right|}$$

(8e)

where $Q_{g_n}$ indicates all successful $g_n$ during each iteration; ${\phi }^{\prime }$ denotes a weighting coefficient.

Because it is possible that the solution $x_2$ does not have a superior objective function to the solution $x_n$ (i.e., $f\left(x_2\right)>f(x_n)$), solution $x_3$ is created to provide the problem a second opportunity to develop an appropriate solution, which is described as follows:

$$\begin{gathered} ifrand < 0.3 \\ {x}_3=x_n+g_n.\left(x_{best}-x_n\right)+0.5.g_n.\left(x_{btr,r1}-x_{por,r2}\right) \\ end \end{gathered}$$

(9)

A dynamic mechanism is used to adjust the population size at each iteration of the EMRUN algorithm. Large quantities of population size ($Np$) used in the early iterations can aid the algorithm in implementing better exploration. In contrast, modest amounts of $Np$ used in the last iteration can assist the algorithm in implementing appropriate exploitation. The dynamic population size that has been suggested can be defined as follows:

$$N{p}^{\prime }=round\left[N{p}_{max}+\left(N{p}_{min}-N{p}_{max}\right).{\left(\frac{Iter}{MaxIter}\right)}^{0.5}\right]$$

(10)

where $N{p}_{max}$ and $N{p}_{min}$ are the maximum and minimum of $Np$, respectively.

This study employs a robust local search mechanism based on the CMA-ES and SQP to boost the best-so-far position at each iteration. This assists the EMRUN in increasing its accuracy and convergence speed. The EMRUN’s pseudocode is expressed in Algorithm 1. According to Algorithm 1, the SQP is applied at the last iterations and with a low chance after $0.2\times Iter$. The EMRUN’s flowchart is depicted in Figure 1.

Figure 1. Flowchart of EMRUN algorithm.

Figure 1. Flowchart of EMRUN algorithm.

Algorithm 1. EMRUN pseudocode

Initialization

Generate the initial population $x_n\left(n=1,2,\dots ,Np\right)$

Calculate the fitness function of the initial population

Determine the solutions $x_{lbest}$ and $x_{best}$ Set $S=0.05$

EMRUN main loop

for $Iter=1: MaxIter$

for $n=1 : Np$

Determine the solutions $x_{btr,r1}$ and $x_{por,r2}$

Calculate $\psi \left(I{d}_n\right)$, $p{c}_{1,n}$, $p{c}_{2,n}$, and $g_n$ using Equations (4d), (5a) and (8a,c).

Calculate the position $x_1$ using Equation (5)

if $f(x_1)<f(x_n)$

$x_{n+1}=x_1$

end

IESQ operator

if $rand<{\left(1-Iter/MaxIter\right)}^2 or Iter<0.1.MaxIter$

Calculate position $x_2$ using Equation (8)

if $f(x_2)<f(x_n)$

$x_{n+1}=x_2$

else

if $rand<0.3$

Calculate position $x_3$ using Equation (9)

if $f(x_3)<f(x_n)$

$x_{n+1}=x_3$

end

end

end

end

if $Iter>0.8.MaxIter or \left(rand\langle S \& Iter\rangle 0.2.MaxIter \right)$

$S=0.1.S$

Apply SQP to $x_{best}$

else

Apply CMA-ES to $x_{best}$

end

Update positions $x_{lbest}$, $x_{bt}$, and $x_{best}$

end for

Update Np using Equation (10)

It = It + 1

end

Phase 3. return $x_{best}$

2.3. Computational Complexity of EMRUN

It is possible to express the computational complexity of an optimization technique as a function that describes the connection between the running time of the algorithm and the size of the issue it is solving. To achieve this objective, the Big-O notation is often used as a primary mode of communication. The components listed below comprise the majority of the suggested EMRUN: initialization, objective function calculation, sorting, and updating the solutions. As a direct result of this analysis, we were able to formulate an estimation of the algorithm’s degree of complexity: 𝛰 [𝑁p × [1 + $MaxIter$× 𝑁p × [1 + 𝑙𝑜𝑔 𝑁p + 𝐷]]].

3. Validation of EMRUN Using Test Functions

As part of this research, the EMRUN algorithm is evaluated with seven metaheuristic algorithms tested against the CEC-2017 test sets, which have 29 test functions with varying features. The seven optimization methods consist of the hybrid adaptive teaching-learning-based optimization (TLBO) with DE (ATLDE) [44], an adaptive DE and PSO (ADEPSO) [45], adaptive differential evolution with optional external archive (JADE) [46], generalized oppositional TLBO (GOTLBO) [47], performance-guided JAYA (PGJAYA) [48], RUN, and Cauchy and gaussian sine cosine optimization (CGSCA) [49].

Table 1. Parameter settings for the eight algorithms.

Algorithms	Parameters
EMRUN	$N=50$, $\omega =0.9$
JADE	$N=50$, $c=0.9$
ATLDE	$N=50$, $CR=0.9$
ADEPSO	$N=50$, $\lambda =0.001$
GOTLBO	$N=50$, Jumping rate $J_r=0.3$
PGJAYA	$N=50$
CGSCA	$N=50$, $a=2$, $S_f=1$, $j=1$
RUN	$N=50$

shows the parameters of these seven optimization methods and EMRUN’s default settings. The maximum number of function evaluations ($MaxNFE$) and the $Np$ is 6000D and 3D, respectively, for this study. The dimensionality of the issue is represented by D, and the functions are reviewed on a scale of 30D. Unimodal functions (UFs), multimodal functions (MFs), hybrid functions (HFs), and composite functions (CFs) are all included in the benchmark test functions (

Table 2. Specifications of the CEC-2017 test functions.

Type	Function No.	Optimum	Feasible Space
UFs	1	100	[−100,100]
	3	300	[−100,100]
MFs	4	400	[−100,100]
	5	500	[−100,100]
	6	600	[−100,100]
	7	700	[−100,100]
	8	800	[−100,100]
	9	900	[−100,100]
	10	1000	[−100,100]
HFs	11	1100	[−100,100]
	12	1200	[−100,100]
	13	1300	[−100,100]
	14	1400	[−100,100]
	15	1500	[−100,100]
	16	1600	[−100,100]
	17	1700	[−100,100]
	18	1800	[−100,100]
	19	1900	[−100,100]
	20	2000	[−100,100]
CFs	21	2100	[−100,100]
	22	2200	[−100,100]
	23	2300	[−100,100]
	24	2400	[−100,100]
	25	2500	[−100,100]
	26	2600	[−100,100]
	27	2700	[−100,100]
	28	2800	[−100,100]
	29	2900	[−100,100]
	30	3000	[−100,100]

).

3.1. Results of the Benchmark Test Functions

The outcomes of the CEC-2017 functions are shown in

Table A1. Mean and SD values achieved by eight algorithms for CEC-17.

Function		EMRUN	RUN	ADEPSO	ATLDE	CGSCA	GOTLBO	JADE	PGJAYA
f1	Mean	1.00 × 102	3.62 × 103	1.00 × 102	1.00 × 102	8.02 × 108	1.82 × 1010	1.00 × 102	3.89 × 103
	SD	1.55 × 10−10	4.91 × 103	8.27 × 10−9	2.85 × 10−2	3.95 × 108	5.19 × 109	1.69 × 10−14	5.42 × 103
f3	Mean	3.00 × 102	3.00 × 102	3.00 × 102	3.80 × 103	1.64 × 103	4.24 × 104	1.68 × 104	3.03 × 102
	SD	7.14 × 10−13	6.52 × 10−2	1.90 × 10−12	2.22 × 103	7.76 × 102	1.25 × 104	2.29 × 104	1.82 × 100
f4	Mean	4.12 × 102	4.59 × 102	4.29 × 102	4.68 × 102	4.53 × 102	3.47 × 103	4.14 × 102	4.19 × 102
	SD	2.34 × 101	3.53 × 101	3.09 × 101	1.54 × 101	2.71 × 101	1.22 × 103	2.43 × 101	2.79 × 101
f5	Mean	5.20 × 102	6.87 × 102	6.56 × 102	5.39 × 102	5.50 × 102	7.29 × 102	5.31 × 102	6.02 × 102
	SD	3.59 × 100	3.79 × 101	3.70 × 101	1.10 × 101	6.08 × 100	2.75 × 101	5.45 × 100	2.08 × 101
f6	Mean	6.00 × 102	6.37 × 102	6.36 × 102	6.00 × 102	6.20 × 102	6.43 × 102	6.00 × 102	6.03 × 102
	SD	9.67 × 10−2	9.54 × 100	7.50 × 100	4.47 × 10−2	4.00 × 100	7.02 × 100	1.83 × 10−7	5.71 × 100
f7	Mean	7.53 × 102	1.01 × 103	1.01 × 103	7.71 × 102	7.78 × 102	1.08 × 103	7.61 × 102	8.49 × 102
	SD	7.60 × 100	6.17 × 101	6.32 × 101	1.11 × 101	9.30 × 100	6.08 × 101	4.46 × 100	3.04 × 101
f8	Mean	8.22 × 102	9.43 × 102	9.26 × 102	8.42 × 102	8.39 × 102	9.94 × 102	8.32 × 102	8.96 × 102
	SD	6.96 × 100	2.30 × 101	2.36 × 101	1.01 × 101	7.92 × 100	2.05 × 101	4.76 × 100	3.05 × 101
f9	Mean	9.01 × 102	3.41 × 103	3.07 × 103	9.04 × 102	1.03 × 103	3.85 × 103	9.00 × 102	2.58 × 103
	SD	2.59 × 100	6.35 × 102	7.43 × 102	1.00 × 101	8.02 × 101	8.33 × 102	1.15 × 10−1	7.17 × 102
f10	Mean	2.69 × 103	4.51 × 103	4.61 × 103	3.55 × 103	2.23 × 103	6.88 × 103	3.33 × 103	4.83 × 103
	SD	2.74 × 102	7.24 × 102	5.97 × 102	1.18 × 103	2.39 × 102	6.02 × 102	2.79 × 102	6.57 × 102
f11	Mean	1.15 × 103	1.19 × 103	1.26 × 103	1.15 × 103	1.20 × 103	2.08 × 103	1.14 × 103	1.18 × 103
	SD	2.55 × 101	3.14 × 101	4.68 × 101	3.45 × 101	3.87 × 101	3.48 × 102	2.94 × 101	3.42 × 101
f12	Mean	2.47 × 103	8.66 × 104	2.04 × 104	3.24 × 104	1.87 × 107	1.35 × 109	3.24 × 103	3.00 × 104
	SD	3.68 × 102	3.59 × 104	1.60 × 104	1.98 × 104	1.26 × 107	9.05 × 108	1.42 × 103	1.34 × 104
f13	Mean	1.37 × 103	2.51 × 104	2.78 × 103	1.46 × 103	4.15 × 104	5.87 × 105	1.34 × 103	1.45 × 104
	SD	4.85 × 101	9.88 × 103	7.61 × 102	2.09 × 102	3.17 × 104	1.82 × 106	2.41 × 101	1.07 × 104
f14	Mean	1.47 × 103	1.89 × 103	1.63 × 103	1.44 × 103	1.69 × 103	1.59 × 103	8.87 × 103	2.36 × 103
	SD	3.03 × 101	3.57 × 102	6.72 × 101	1.43 × 101	2.93 × 102	4.69 × 101	1.27 × 104	9.25 × 102
f15	Mean	1.56 × 103	1.41 × 104	1.78 × 103	1.54 × 103	2.32 × 103	1.35 × 104	2.51 × 103	6.99 × 103
	SD	2.83 × 101	2.91 × 103	1.29 × 102	3.59 × 101	5.66 × 102	8.82 × 103	3.40 × 103	5.09 × 103
f16	Mean	1.89 × 103	2.58 × 103	2.57 × 103	2.12 × 103	1.73 × 103	2.86 × 103	2.18 × 103	2.43 × 103
	SD	1.80 × 102	2.85 × 102	2.28 × 102	2.03 × 102	6.33 × 101	2.68 × 102	1.43 × 102	2.37 × 102
f17	Mean	1.75 × 103	2.12 × 103	2.30 × 103	1.83 × 103	1.78 × 103	2.12 × 103	1.80 × 103	2.05 × 103
	SD	3.08 × 101	2.04 × 102	1.57 × 102	1.01 × 102	1.48 × 101	1.56 × 102	3.51 × 101	2.00 × 102
f18	Mean	1.88 × 103	2.38 × 104	7.37 × 103	2.06 × 103	2.12 × 105	4.08 × 104	3.03 × 104	5.42 × 104
	SD	3.25 × 101	1.13 × 104	9.14 × 103	4.08 × 102	1.93 × 105	2.87 × 104	5.44 × 104	2.65 × 104
f19	Mean	1.95 × 103	6.30 × 103	2.05 × 103	1.92 × 103	3.52 × 103	9.43 × 104	2.90 × 103	9.30 × 103
	SD	2.39 × 101	1.48 × 103	6.58 × 101	2.11 × 101	2.27 × 103	1.34 × 105	3.92 × 103	7.24 × 103
f20	Mean	2.08 × 103	2.39 × 103	2.45 × 103	2.14 × 103	2.10 × 103	2.47 × 103	2.15 × 103	2.33 × 103
	SD	5.80 × 101	1.34 × 102	1.63 × 102	1.13 × 102	2.20 × 101	1.32 × 102	5.61 × 101	1.78 × 102
f21	Mean	2.32 × 103	2.43 × 103	2.44 × 103	2.34 × 103	2.22 × 103	2.51 × 103	2.33 × 103	2.40 × 103
	SD	4.87 × 100	3.79 × 101	2.84 × 101	8.39 × 100	2.65 × 101	2.48 × 101	6.02 × 100	2.52 × 101
f22	Mean	2.30 × 103	3.05 × 103	2.40 × 103	3.02 × 103	2.37 × 103	3.88 × 103	2.30 × 103	2.97 × 103
	SD	6.00 × 10−1	1.56 × 103	5.68 × 102	1.28 × 103	3.08 × 101	5.85 × 102	0.00 × 100	1.54 × 103
f23	Mean	2.67 × 103	2.77 × 103	2.93 × 103	2.69 × 103	2.65 × 103	3.00 × 103	2.67 × 103	2.77 × 103
	SD	7.88 × 100	3.43 × 101	9.38 × 101	5.54 × 101	8.15 × 100	7.14 × 101	5.95 × 100	3.10 × 101
f24	Mean	2.84 × 103	2.91 × 103	3.10 × 103	2.87 × 103	2.72 × 103	3.28 × 103	2.84 × 103	2.92 × 103
	SD	1.03 × 101	2.69 × 101	1.05 × 102	1.19 × 101	9.10 × 101	8.33 × 101	6.01 × 100	3.43 × 101
f25	Mean	2.89 × 103	2.90 × 103	2.89 × 103	2.89 × 103	2.96 × 103	3.51 × 103	2.89 × 103	2.89 × 103
	SD	2.38 × 100	1.36 × 101	9.84 × 100	1.18 × 100	1.42 × 101	1.91 × 102	6.98 × 10−1	1.02 × 101
f26	Mean	3.70 × 103	4.81 × 103	5.63 × 103	4.09 × 103	3.09 × 103	6.89 × 103	3.77 × 103	4.11 × 103
	SD	3.26 × 102	1.23 × 103	1.54 × 103	1.80 × 102	5.33 × 101	9.18 × 102	1.76 × 102	1.09 × 103
f27	Mean	3.22 × 103	3.26 × 103	3.33 × 103	3.22 × 103	3.10 × 103	3.51 × 103	3.21 × 103	3.23 × 103
	SD	1.20 × 101	1.84 × 101	9.92 × 101	1.56 × 101	2.07 × 100	1.05 × 102	4.76 × 100	1.66 × 101
f28	Mean	3.12 × 103	3.11 × 103	3.14 × 103	3.17 × 103	3.27 × 103	4.30 × 103	3.14 × 103	3.14 × 103
	SD	4.69 × 101	4.03 × 101	6.68 × 101	6.83 × 101	5.24 × 101	3.29 × 102	5.80 × 101	5.34 × 101
f29	Mean	3.36 × 103	3.96 × 103	4.10 × 103	3.46 × 103	3.24 × 103	4.43 × 103	3.39 × 103	3.58 × 103
	SD	3.00 × 101	2.58 × 102	3.49 × 102	1.26 × 102	3.74 × 101	2.16 × 102	2.02 × 101	1.61 × 102
f30	Mean	5.57 × 103	2.97 × 104	5.83 × 103	5.22 × 103	8.04 × 105	3.28 × 106	5.26 × 103	8.42 × 103
	SD	8.70 × 102	1.46 × 104	1.14 × 103	2.57 × 102	4.09 × 105	3.01 × 106	1.83 × 102	3.28 × 103

(see Appendix A). It can be noticed that EMRUN surpasses RUN, ADEPSO, JADE, CGSCA, ATLDE, PGJAYA, and GOTLBO on functions 3, 4, 5, 7, 8, 12, 17, 18, 20, 22, 26, and 29 of the CEC 2017 set, proving that the EMRUN algorithm suggested in this study beats all seven other methods. The following is a summary of the relevant details:

Unimodal functions (UFs: f₁ and f₃): Table A1 shows that EMRUN, ATLDE, ADEPSO, and JADE have the most excellent performance on function 1 in terms of mean, while EMRUN has the best efficiency on function 3 in terms of mean (3.00 × 10²) and SD (7.14 × 10⁻¹³). The GOTLBO has the worst performance on functions 1 and 3. According to the findings, the CMA-ES, SQP, and IESQ can easily capture fitness landscape features in the UFs and then assist a population in seeking attractive locations.

Multimodal functions (MFs: f₄ and f₁₀): According to the comparison findings on the MFs shown in Table A1, EMRUN achieves the best result on functions 4, 5, 6, 7, and 8, while JADE achieves the best result on functions 6 and 9. In addition, CGSCA has the best efficiency on function 10. Even though EMRUN, along with the other seven rivals, cannot obtain the best result on test functions 9 and 10, it achieves the second-best rank on these functions. As shown by the stable and promising feature of EMRUN, it is possible to achieve a balance between maintaining population variety and speeding up convergence on certain MFs using the suggested RK-SM operator, CMA-ES, and adaptive parameters.

Hybrid functions (HFs: f₁₁–f₂₀) and composite functions (CFs: f₂₁–f₃₀): HFs and CFs are two forms of sophisticated functions that are exceedingly challenging. The findings on f11–f30 reported in Table A1 demonstrate that the algorithms MERUN, ATLDE, and JADE outperform their counterparts in terms of competitiveness. In the case of HFs, EMRUN achieves the most outstanding results on functions 12, 17, 18, and 20. Furthermore, the suggested technique produces the best results for the CFs on functions 22, 26, and 29. Although EMRUN cannot obtain top results on other test functions, it ranks second. The comparative findings show that the adaptive parameters RK-SM, IESQ, and CMA-ES are useful for balancing exploration and exploitation, allowing the population to have stable search behavior.

3.2. Ranking Analysis of EMRUN

The Wilcoxon rank (WR) test [50] was used to examine whether the differences between the results supplied by the EMRUN and a particular algorithm are statistically crucial at the 5% level of significance, which is included in

Table A2. Results of the WT achieved by eight algorithms on CEC-17.

F	EMRUN vs. ADEPSO				EMRUN vs. ATLDE				EMRUN vs. CGSCA				EMRUN vs. GOTLBO
	S+	S−	P−V	W	S+	S−	P−V	W	S+	S−	P−V	W	S+	S−	P−V	W
1	439	26	2.34 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
3	461.5	3.5	2.65 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
4	374	91	3.89 × 10−3	+	463	2	2.13 × 10−6	+	445	20	1.24 × 10−5	+	465	0	1.73 × 10−6	+
5	465	0	1.73 × 10−6	+	458	7	3.52 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
6	465	0	1.73 × 10−6	+	218	247	7.66 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
7	465	0	1.73 × 10−6	+	464	1	1.92 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	−
8	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+	461	4	2.60 × 10−6	+	465	0	1.73 × 10−6	+
9	465	0	1.73 × 10−6	+	311	154	1.04 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
10	465	0	1.73 × 10−6	+	406	59	3.59 × 10−4	+	11	454	5.22 × 10−6	−	465	0	1.73 × 10−6	+
11	465	0	1.73 × 10−6	+	248	217	7.50 × 10−1	=	450	15	7.69 × 10−6	+	465	0	1.73 × 10−6	+
12	465	0	1.73 × 10−6	+	464	1	1.92 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
13	465	0	1.73 × 10−6	+	290	175	2.37 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
14	465	0	1.73 × 10−6	+	17	448	9.32 × 10−6	−	456	9	4.29 × 10−6	+	465	0	1.73 × 10−6	+
15	465	0	1.73 × 10−6	+	160	305	1.36 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
16	465	0	1.73 × 10−6	+	431	34	4.45 × 10−5	+	83	382	2.11 × 10−3	−	465	0	1.73 × 10−6	+
17	465	0	1.73 × 10−6	+	400	65	5.71 × 10−4	+	401	64	5.29 × 10−4	+	465	0	1.73 × 10−6	+
18	465	0	1.73 × 10−6	+	444	21	1.36 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
19	455	10	4.73 × 10−6	+	60	405	3.88 × 10−4	−	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
20	465	0	1.73 × 10−6	+	356	109	1.11 × 10−2	+	363	102	7.27 × 10−3	+	465	0	1.73 × 10−6	+
21	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+	1	464	1.92 × 10−6	−	465	0	1.73 × 10−6	+
22	349.5	115.5	2.67 × 10−2	+	384	81	1.83 × 10−3	+	464	1	1.92 × 10−6	+	465	0	1.73 × 10−6	+
23	465	0	1.73 × 10−6	+	435	30	3.11 × 10−5	+	1	464	1.92 × 10−6	−	465	0	1.73 × 10−6	+
24	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+	0	465	1.73 × 10−6	−	465	0	1.73 × 10−6	+
25	273	192	4.05 × 10−1	=	183	282	3.09 × 10−1	=	464	1	1.92 × 10−6	+	465	0	1.73 × 10−6	+
26	447	18	1.02 × 10−5	+	446	19	1.13 × 10−5	+	10	455	4.73 × 10−6	−	465	0	1.73 × 10−6	+
27	465	0	1.73 × 10−6	+	303	162	1.47 × 10−1	=	0	465	1.73 × 10−6	−	465	0	1.73 × 10−6	+
28	326	139	4.93 × 10−2	+	436	29	2.84 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
29	465	0	1.73 × 10−6	+	411	54	2.41 × 10−4	+	0	465	1.73 × 10−6	−	465	0	1.73 × 10−6	+
30	331	134	4.28 × 10−2	+	131	334	3.68 × 10−2	−	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
F	EMRUN vs. JADE				EMRUN vs. PGJAYA				EMRUN vs. RUN
	S+	S−	P−V	W	S+	S−	P−V	W	S+	S−	P−V	W
1	0	465	1.45 × 10−7	−	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
3	275	190	1.95 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
4	244	221	8.13 × 10−1	=	332	133	4.07 × 10−2	+	453	12	5.75 × 10−6	+
5	431	34	4.45 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
6	0	465	1.73 × 10−6	−	378	87	2.77 × 10−3	+	465	0	1.73 × 10−6	+
7	437	28	2.60 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
8	444	21	1.36 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
9	9	456	4.25 × 10−6	−	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
10	447	18	1.02 × 10−5	+	465	0	1.73 × 10−6	+	464	1	1.92 × 10−6	+
11	196	269	4.53 × 10−1	=	412	53	2.22 × 10−4	+	453	12	5.75 × 10−6	+
12	329	136	4.72 × 10−2	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
13	66	399	6.16 × 10−4	−	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
14	255	210	6.44 × 10−1	=	463	2	2.13 × 10−6	+	465	0	1.73 × 10−6	+
15	140	325	5.71 × 10−2	=	456	9	4.29 × 10−6	+	465	0	1.73 × 10−6	+
16	458	7	3.52 × 10−6	+	460	5	2.88 × 10−6	+	465	0	1.73 × 10−6	+
17	451	14	6.98 × 10−6	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
18	236	229	9.43 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
19	64	401	5.29 × 10−4	−	463	2	2.13 × 10−6	+	465	0	1.73 × 10−6	+
20	424	41	8.19 × 10−5	+	464	1	1.92 × 10−6	+	465	0	1.73 × 10−6	+
21	445	20	1.24 × 10−5	+	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
22	0	465	2.28 × 10−5	−	351	114	1.48 × 10−2	+	399	66	6.16 × 10−4	+
23	303	162	1.47 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
24	229	236	9.43 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+
25	159	306	1.31 × 10−1	=	84	381	2.26 × 10−3	−	393	72	9.63 × 10−4	+
26	214	251	7.04 × 10−1	=	315	150	8.97 × 10−2	=	414	51	1.89 × 10−4	+
27	76	389	1.29 × 10−3	−	454	11	5.22 × 10−6	+	465	0	1.73 × 10−6	+
28	213	252	7.12 × 10−1	=	416	49	1.60 × 10−4	+	359	106	9.27 × 10−3	+
29	334	131	3.68 × 10−2	+	453	12	5.75 × 10−6	+	465	0	1.73 × 10−6	+
30	162	303	1.47 × 10−1	=	465	0	1.73 × 10−6	+	465	0	1.73 × 10−6	+

(see Appendix A). Each optimization method’s results are shown with a ‘+’, ‘−’ or ‘=’, meaning that it is worse than, better than, or equal to the EMRUN. According to the WR test (

Table 3. Summary results of WT achieved by EMRUN vs. other methods.

EMRUN vs.	+/−/=
ADEPSO	28/0/1
ATLDE	19/3/7
CGSCA	21/8/0
GOTLBO	28/1/0
JADE	10/7/12
PGJAYA	27/1/1
RUN	29/0/0

), EMRUN performs significantly better than ADEPSO, ATLDE, CGSCA, GOTLBO, JADE, PGJAYA, and RUN on 28, 19, 21, 28, 10, 27, and 29 test functions and somewhat worse on 0, 3, 8, 1, 7, 1, and 0 test functions. These findings revealed that the EMRUN approach was capable of delivering extraordinary outcomes.

3.3. Convergence Analysis of EMRUN

In this section, to examine the convergence behavior of all optimization methods, a set of 14 test functions is employed. Figure 2 displays the convergence curves computed as the mean of the objective function found in 30 various runs of all methods. Three distinct stages can be seen in the EMRUN curves: exploration (which takes up one-third of the optimization procedure), exploration with exploitation (which takes up the remaining one-third of the optimization procedure), and, finally, exploitation (which takes place after the last two-thirds). According to these phases, two behaviors can be seen in this figure, obtained from the optimization operators such as the RK-SM, IESQ, adaptive parameters, CMA-ES, and SQP. When observing the initial behavior of the suggested EMRUN algorithm, it proceeds toward the optimum point with a quick convergence rate via the IESQ, the CMA-ES, and the adaptive parameters, indicating adequate performance for optimizing issues in this context. This behavior is also evident in functions 1, 3, 4, 9, 12, 14, 19, 22, and 30. In the second case, the EMRUN begins the optimization process with a sluggish convergence speed and progressively speeds up as it nears the ideal point. The RK-SM, CMA-ES, and IESQ are all responsible for this behavior. Further evidence of this tendency is seen in functions 5, 7, 17, 25, and 28. As a result, compared with competitive algorithms, the suggested algorithm, which uses the operators of RK-SM, IESQ, CMA-ES, SQP, and adaptive parameters, has a more robust capacity to leap out of local optima when it finds itself in one of them.

4. Modeling the Nonlinear Muskingum Flood Routing

Using the continuity and the storage–discharge relationship, the generalized nonlinear Muskingum (GNLM4) technique of flood routing introduced by [51] can be implemented, which is formulated as follows:

$$\frac{d{s}_t}{dt}\approx \frac{\Delta s_t}{\Delta t}=I_t-O_t$$

(11)

where $s_t$ expresses the storage volume, $I_t$ expresses the inflow, $O_t$ indicates outflow in a river, and $t$ expresses the time step. The relationship of GNLM4 can be formulated as follows:

$$s_t=K\left[X_1\left(C_1{I}_t^{m1}\right)+X_2\left(C_1{I}_{t+1}^{m1}\right)+\left(1-X_1-X_2\right)\left(C_2{O}_t^{m2}\right)\right]{ }^{\alpha }$$

(12a)

where $K$ expresses a time-storage coefficient, $X_1$ and $X_2$ denote the weighting coefficients for determining the importance of the inflow in the $t^{th}$ and ${\left(t+1\right)}^{th}$ time steps, respectively. $C_1$ and $C_2$ indicate the factors of inflow and outflow, respectively.

The main steps of the GNLM4 model are expressed as follows:

(1) Estimate the main parameters of GNLM4 ($K$, $X_1$, $X_2$, $C_1$, $C_2$, $m1$, $m2$, and $\alpha$).

(2) Compute the initial storage $s_0$ using Equation (12) under the assumption that the initial outflow is equal to the initial inflow (${\widehat{O}}_0=I_0$):

$$s_0=K\left[X_1\left(C_1{I}_0^{m1}\right)+X_2\left(C_1{I}_1^{m1}\right)+\left(1-X_1-X_2\right)\left(C_2{O}_0^{m2}\right)\right]{ }^{\alpha }$$

(12b)

(3) Compute the rate of variation in the $s_t$ at the $t^{th}$ time step, ($\frac{\Delta s_t}{\Delta t}$ ) with Equation (13):

$$\frac{\Delta s_t}{\Delta t}=I_t-{\left(\frac{1}{C_2}\left(\left[\frac{1}{1-X_1-X_2}\right].{\left(\frac{s_t}{K}\right)}^{\frac{1}{\alpha }}-\left[\frac{X_1}{1-X_1-X_2}\right].C_1{I}_t^{m1}-\left[\frac{X_2}{1-X_1-X_2}\right].C_1{I}_{t+1}^{m1}\right)\right)}^{1/m2}$$

(13)

where $m1$, $m2$, and $\alpha$ denote parameters for specifying the flood nonlinearity.

(4) Compute the $s_{t+1}$, with Equation (14):

$$s_{t+1}=s_t+\Delta s_t$$

(14)

(5) Compute the ${\widehat{O}}_t$ with Equation (15):

$${\widehat{O}}_t={\left(\frac{1}{C_2}\left(\left[\frac{1}{1-X_1-X_2}\right].{\left(\frac{s_t}{K}\right)}^{\frac{1}{\alpha }}-\left[\frac{X_1}{1-X_1-X_2}\right].C_1{I}_t^{m1}-\left[\frac{X_2}{1-X_1-X_2}\right].C_1{I}_{t+1}^{m1}\right)\right)}^{1/m2}$$

(15)

(6) Steps (3) through (5) should be repeated until the simulation has been completed.

The fitness function for estimating the parameters of the GNLM4 model was set as Equation (16):

$$Minimize SSA={\displaystyle \sum }_{t=1}^T{\left({\widehat{O}}_t-O_t\right)}^2$$

(16)

where $SSA$ denotes the sum of squared deviations.

Furthermore, the results calculated by the optimization methods are assessed by various error metrics comprising the deviations of the peak of outflows (DPO), the mean absolute relative error (MARE), the absolute deviations (SAD), and the variance explained in a percentage ($VarexQ$) [1]. These error metrics are defined as follows:

$$DPO=\left|{\widehat{O}}_t{}^{peak}-O_t{}^{peak}\right|$$

(17)

$$MARE=\frac{1}{T}{\displaystyle \sum }_{t=1}^T\frac{\left|{\widehat{O}}_t-O_t\right|}{{\widehat{O}}_t}$$

(18)

$$SAD={\displaystyle \sum }_{t=1}^T\left|{\widehat{O}}_t-O_t\right|$$

(19)

$$VarexQ=\left[1-\frac{{{\displaystyle \sum }}_{t=1}^T{\left({\widehat{O}}_t-O_t\right)}^2}{{{\displaystyle \sum }}_{t=1}^T{\left(O_t-O^{avg}\right)}^2}\right]\times 100$$

(20)

where ${\widehat{O}}_t{}^{peak}$ and $O_t^{peak}$ express the peak flow amounts of the observed and computed outflows and $O_t^{avg}$ expresses the average of the outflows.

4.1. Results of the GNLM4 Models

The proposed EMRUN is employed in this study to optimize the GNLM4 parameters. The problems have a nonlinear and nonconvex structure; thus, they are challenging enough for any optimization method. The proposed method is applied to two case studies (a smooth single-peak hydrograph and a multimodal example). The following sections compare the proposed EMRUN with seven optimization methods in all case studies.

4.2. First Case Study: Wilson Flood (Single-Peak Hydrograph)

As a first case study, the Wilson flood’s hydrographs (WF) [52] (i.e., the inflow and outflow hydrographs) were used, where the number of time steps ($T$) and the length of the simulation time step ($\Delta t$) are equal to 21 and 6, respectively. The bounds of the unknown parameters to optimize this problem are $0.1\le K\le 10$, $0\le X_1\&X_2\le 1$, $0.1\le C_1\&C_2\le 2$, $1\le \alpha \le 5$, and $0.1\le m1\&m2\le 1$.

Table 4. Best, AVG, worst, and SD values achieved by eight algorithms on the WF problem.

	Optimization Methods
	ADEPSO	ATLDE	CGSCA	GOTLBO	EMRUN	JADE	PGJAYA	RUN
Best	4.76	4.91	6.43	6.17	4.76	4.76	17.58	4.78
AVG	5.07	6.42	25.55	15.14	4.76	13.82	32.25	5.48
Worst	6.70	10.27	57.30	50.01	4.76	32.04	42.33	7.64
SD	6.35 × 10−1	1.50 × 100	1.55 × 101	1.06 × 101	5.87 × 106	1.05 × 101	6.67 × 100	6.31 × 101

reports the statistical results of all optimization methods in terms of the best, average (AVG), worst, and standard deviation (SD) amounts of fitness function over 30 various runs. From the table, the EMRUN, ADEPSO, and JADE can achieve better performance in terms of the best value of the objective function. In addition, the EMRUN algorithm can outperform all contestant optimization methods in terms of the worst (4.76), AVG (4.76), and SD (5.87 × 10⁶) amounts.

Table 5. Eight algorithms obtained four error metrics on the WF problem.

	Optimization Methods
	ADEPSO	ATLDE	CGSCA	GOTLBO	EMRUN	JADE	PGJAYA	RUN
SAD	6.720	6.593	8.496	8.579	6.654	6.654	15.555	6.715
DPO	0.090	0.014	0.112	0.190	0.081	0.081	0.488	0.114
MARE	0.011	0.011	0.013	0.014	0.011	0.011	0.021	0.011
VarexQ	99.961	99.960	99.947	99.949	99.961	99.961	99.856	99.961

lists four statistical metrics (i.e., $DPO$, $MARE$, $SAD$, and $VarexQ$) to further assess the efficiency of the proposed algorithm compared with the other methods. According to the table, the EMRUM and JADE algorithm achieved the best amounts for the $DPO$, $MARE$, $SAD$, and $VarexQ$ when compared with the other methods. These superior results come from the effective operators and mechanisms implemented in the proposed EMRUN algorithm, including the RK-SM and the IESQ, CMA-ES, and SQP methods. Figure 3 displays the convergence graphs of all methods for the problem. This figure is based on the average fitness function over 30 various runs. The figure shows that the EMRUM algorithm has a superior convergence speed to the other algorithms, followed by ADEPSO, ATLDE, RUN, GOTLBO, JADE, CGSCA, and PGJAYA, in order.

4.3. Second Case Study: Viessman and Lewis Flood (Multipeak Hydrograph)

As a second case study, Viessman and Lewis’s hydrographs [53] (VLF) were used, where the number of time steps ($T$) and the length of the simulation time step ($\Delta t$) are equal to 23 and 1, respectively. The bounds of the unknown parameters to optimize this problem are $0.1\le K\le 10$, $-1\le X_1\&X_2\le 1$, $0.1\le C_1\&C_2\le 2$, $1\le \alpha \le 5$, and $0.1\le m1\&m2\le 2$. The statistical results of all optimization techniques are shown in

Table 6. Best, AVG, worst, and SD values achieved by eight algorithms on the VLF problem.

	Optimization Methods
	ADEPSO	ATLDE	CGSCA	GOTLBO	EMRUN	JADE	PGJAYA	RUN
Best	52,456.99	53,147.89	71,124.74	56,829.88	52,456.84	56,330.10	66,882.10	53,477.51
AVG	59,640.19	64,453.07	84,879.90	67,071.71	52,456.84	62,176.57	79,825.56	77,492.05
Worst	82,227.76	82,223.41	89,203.02	82,469.00	52,456.89	67,710.32	90,227.46	82,223.65
SD	11,305.195	10,798.381	4331.271	7203.845	0.011	3430.916	6219.807	9519.628

in terms of the best, AVG, worst, and SD amounts of fitness function for 30 different runs. According to the table, the EMRUN may obtain the most excellent performance in terms of the best (52,456.84), worst (52,456.89), AVG (52,456.84), and SD (0.011) quantities of objective function during 30 various runs.

Table 7. Eight algorithms obtained four error metrics on the VLF problem.

	Optimization Methods
	ADEPSO	ATLDE	CGSCA	GOTLBO	EMRUN	JADE	PGJAYA	RUN
SAD	890.26	911.01	1024.03	943.84	889.41	934.89	1010.89	914.18
DPO	9.80	18.79	81.41	23.90	9.66	18.33	55.77	15.52
MARE	0.08	0.08	0.08	0.08	0.08	0.08	0.08	0.08
VarexQ	98.793	98.777	98.364	98.692	98.793	98.704	98.461	98.770

contains four statistical measures (i.e., $DPO$, $MARE$, $SAD$, and $VarexQ$) that can be used to assess the effectiveness of the EMRUN in comparison with the other approaches. According to the table, the EMRUM algorithm outperformed the other approaches in terms of $DPO$ (9.66), $MARE$ (0.08), $SAD$ (889.41), and $VarexQ$ (98.793). These excellent results are due to the efficient operators and mechanisms used in the proposed EMRUN algorithm, including the RK-SM, IESQ, adaptive parameters and the SQP method. An illustration of all the approaches’ convergent graphs is shown in Figure 4. It is calculated using an average fitness function across a total of 30 different runs. According to the graphic, the EMRUM algorithm has a faster convergence speed than the other algorithms, followed by the ADEPSO, JADE, ATLDE, GOTLBO, RUN, PGJAYA, and CGSCA algorithms, in order.

5. Multireservoir Problem

An evaluation of the efficiency and the capabilities of EMRUN was carried out in this area using a six-reservoir system. An improved version of the widely tested multireservoir problem introduced by [54] is presented here. A schematic representation of this system is shown in Figure 5, where the system contains six reservoirs that generate hydroelectric power in 12 time increments. In reality, the issue in this research was created to become a real-world problem. In addition, the figure shows that there are two distinct reservoir configurations: parallel and series. The inflow to reservoirs for the problem is similar to those introduced by [54]. The fitness function and all constraints of the problem are formulated as follows:

$$Max F={\displaystyle \sum }_{l=1}^L{\displaystyle \sum }_{i=1}^T\left(\frac{PO{W}_i^l}{PO{W}_{max}^l}\right)$$

(21)

$$W_{i+1}^l=W_i^l+I_i^l-O_i^l$$

(22)

$$PO{W}_i^l=\min \left[\frac{g.\delta .{\epsilon }^l.O_i^l.{\overline{H}}_i^l}{1000},PO{W}_{max}^l\right]$$

(23)

$${\overline{H}}_i^l=\frac{h_i^l+h_{i+1}^l}{2}$$

(24)

$$h_i^l=a+b.W_i^l+c.{\left(W_i^l\right)}^2+d.{\left(W_i^l\right)}^3$$

(25)

$$W_{min}^l\le W_i^l\le W_{max}^l$$

(26)

$$O_{min}^l\le O_i^l\le O_{max}^l$$

(27)

where L expresses the total number of reservoirs; T expresses the total number of time steps; $P{W}_{max}^k$ indicates the install capacity (IC) of the $l^{th}$ reservoir; $PO{W}_i^l$ indicates the power produced from the $l^{th}$ reservoir; $g$ is equal to 9.81; $\delta$ denotes a positive number to convert the inflow into volume; ${\epsilon }^l$ is the efficiency of the $l^{th}$ reservoir and is equal to 0.9 for all reservoirs; $W_i^l$ denotes the volume of storage for the $l^{th}$ reservoir; $I_i^l$ and $O_i^l$ indicate the volume of inflow and outflow from the $l^{th}$ reservoir; ${\overline{H}}_i^l$ indicates the average level of water in the $l^{th}$ reservoir; $h_i^l$ expresses the water level in the $l^{th}$ reservoir; $W_{min}^l$ and $W_{max}^l$ indicate the minimum and maximum volumes of storage for the $l^{th}$ reservoir, respectively; and $O_{max}^l$ and denote the highest and lowest amounts of released water from the $l^{th}$ reservoir, respectively. The amounts of c, d, and $PO{W}_i^l$ are reported in

Table 8. Values of parameters for Equation (25).

Coefficients	Reservoirs
	1	2	3	4	5	6
$c$	0.1	0.1	0.1	0.1	−0.3	−0.3
$d$	−0.02	−0.02	−0.02	−0.02	0.005	0.005
$PO{W}_i^l$	1000	1000	1000	1000	5000	5000

for all reservoirs. In addition, the values of parameters a and b are equal to 100 and 10, respectively, for all reservoirs. The problem is nonlinear with several decision variables, so it is very challenging for any optimization method.

The reservoir releases ($O_i^l$) are considered as the decision parameters, in this study. The penalty function (PNFs) in this problem is derived from the limitations specified for the quantities of reservoir storage. The fitness function provided in Equation (21) is rewritten as follows after applying the PNFs:

$$PN{F}_1=\{\begin{array}{c}{\omega }_1.{[W_{T+1}^l-W_1^l]}^2 if V_{T+1}^n\ne V_1^n \\ 0 other wise \end{array}$$

(28)

$$PN{F}_2=\{\begin{array}{c}{\omega }_2.[1-\frac{W_i^l}{W_{min}^l}] if W_i^l<W_{min}^l \\ {\omega }_3.[1-\frac{W_{max}^l}{W_i^l}] if W_i^l>W_{max}^l \\ 0 otherwise \end{array}$$

(29)

where ${\omega }_1$, ${\omega }_2$, and ${\omega }_3$ indicate the coefficients of PNFs and are equal to 20. Consequently, the fitness function is rewritten as follows:

$$Max F={\displaystyle \sum }_{l=1}^L{\displaystyle \sum }_{i=1}^T\left(\frac{PO{W}_i^l}{PO{W}_{max}^l}\right)-PN{F}_1+PN{F}_2)$$

(30)

5.1. Results of Six-Reservoir Problem

The six-reservoir problem was solved by the LINGO 8.0 software to determine the optimum solution, and the global solution was equal to 25.514. In addition, the problem was solved by eight optimization methods: EMRUN, RUN, ADEPSO, ATLDE, CGSCA, GOTLBO, JADE, and PGJAYA. The population size ($Np$) and the total number of function evaluations ($MaxNFE$) were equal to 100 and 180,000, respectively.

Table 9. Best, AVG, worst, and SD values achieved by eight algorithms on the six-reservoir problem.

	Optimization Methods
	ADEPSO	ATLDE	GOTLBO	RUN	JADE	PGJAYA	EMRUN	CGSCA
Best	24.62	22.12	23.08	24.08	24.96	22.80	25.50	22.37
AVG	24.31	21.14	22.52	23.16	24.67	22.04	25.46	18.47
Worst	23.43	19.19	22.10	21.28	24.34	21.46	25.34	9.99
SD	0.215	0.712	0.245	0.659	0.136	0.340	0.042	3.244

lists the statistical outcomes achieved by all algorithms for the six-reservoir problem, comprising the best, worst, AVG, and SD values over 30 various runs. EMRUN’s best, AVG, worst, and SD values are 25.50, 25.46, 25.34, and 0.042, respectively, suggesting that EMRUN performed much better at the task than the other optimizers. The optimum solution value is 25.514; however, the best value from EMRUN is 25.505, which is better than the other methods. Therefore, the EMRUN has the best performance in the problem, followed by the JADE, ADEPSO, RUN, GOTLBO, PGJAYA, ATLDE, and CGSCA methods, in order.

The box plots of the fitness functions for all optimization methods are shown in Figure 6. Patterns in box plots are useful for evaluating distributions of a large number of groups or datasets. Box plots represent the distributions of the objective function throughout the 30 runs for each optimization strategy. From Figure 6, it can be seen that the maximum and minimum amounts of the EMRUN box plot are relatively near to each other and that the deviation of the data is significantly smaller than that of the other techniques. According to the distribution achieved for each algorithm over 30 runs, the EMRUN is more reliable and precise than the other methods. It should be noted that the JADE and ADEPSO have the second-and third-best performance, respectively, in solving the problem.

EMRUN’s convergence graphs are shown in Figure 7, revealing that its precision and speed are superior to its rivals. According to the figure, EMRUN is finding the best solution while other optimization methods are just beginning to explore some solutions. Furthermore, during the early iterations, EMRUN rapidly discovers the best solution, while some of the examined algorithms have slipped into a local optimum. These findings suggest that the EMRUN has a significant capacity to explore and has exploration tendencies and that the two phases exhibit outstanding stability thanks to the effects of the $\psi$ parameter and adaptive parameters (i.e., $p{c}_n$ and $g_n$). Consequently, it can be stated that the EMRUN algorithm can produce accurate and trustworthy results in a shorter amount of time than the other algorithms.

To assess the efficiency of eight optimization methods, eight statistical metrics are used in this work, comprising the correlation coefficient (R), mean absolute percentage error (MAPE), root mean square error (RMSE), Willmott’s agreement index (I_A), median absolute error (MdAE), maximum absolute error (MaxAE), wave hedges distance (WHD), Nash–Sutcliffe model efficiency coefficient (NSE) [55,56], which are formulated, in order, as follows:

$$R=\frac{{{\displaystyle \sum }}_{i=1}^T\left(O_{LG,i}-\overline{O_{LG}}\right).\left(O_{AL,i}-\overline{O_{AL}}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^T{\left(O_{LG,i}-\overline{O_{LG}}\right)}^2{{\displaystyle \sum }}_{i=1}^T{\left(O_{AL,i}-\overline{O_{AL}}\right)}^2}}$$

(31)

$$RMSE={\left(\frac{1}{T}{\displaystyle \sum }_{i=1}^T{\left(O_{LG,i}-O_{AL,i}\right)}^2\right)}^{0.5}$$

(32)

$$MAPE(\%)=\left(\frac{100}{T}\right){\displaystyle \sum }_{i=1}^T\left|\frac{O_{LG,i}-O_{AL,i}}{O_{LG,i}}\right|$$

(33)

$$I_A=1-\frac{{{\displaystyle \sum }}_i^T{\left(O_{LG,i}-O_{AL,i}\right)}^2}{{{\displaystyle \sum }}_{i=1}^T{\left(\left|\left(O_{LG,i}-\overline{O_{LG}}\right)\right|+\left|\left(O_{AL,i}-\overline{O_{LG}}\right)\right|\right)}^2}, 0<I_A\le 1$$

(34)

$$MdAE=media{n}_{i=1,\dots ,T}\left|O_{LG,i}-O_{AL,i}\right|$$

(35)

$$MaxAE=ma{x}_{i=1,\dots ,T}\left|O_{LG,i}-O_{AL,i}\right|$$

(36)

$$WHD={\displaystyle \sum }_{i=1}^T\left(\frac{|O_{LG,i}-O_{AL,i}|}{ma{x}_{i=1,\dots ,T}\left(O_{LG,i},O_{AL,i}\right)}\right)$$

(37)

$$NSE=1-\frac{{{\displaystyle \sum }}_i^T{\left(O_{LG,i}-O_{AL,i}\right)}^2}{{{\displaystyle \sum }}_i^T{\left(O_{LG,i}-{\overline{O}}_{LG}\right)}^2}$$

(38)

where $O_{LG,i}$ expresses the release calculated by LINGO; $O_{AL,i}$ expresses the release calculated by the algorithms; and $\overline{O_{AL}}$ and $\overline{O_{LG}}$ are the average amounts of release, respectively, from LINGO the algorithms.

Three indices consisting of the root mean square error (cRMSE), SD, and R are utilized to obtain the Taylor diagram. The geometric distance and the reference point (target) in the polar area are used to illustrate the ability of the methods [57]. The $cRMSE$ is expressed as follows:

$$cRMS{E}^2=S{D}_{LG}^2+S{D}_{AL}^2-2.R.S{D}_{LG}.S{D}_{AL}$$

(39)

where $S{D}_{LG}$ and $S{D}_{AL}$ indicate the standard deviation of the releases calculated by the LINGO and the algorithms, respectively.

Table 10. Eight algorithms obtained eight error metrics on the six-reservoir problem.

	Optimization Methods
	ADEPSO	ATLDE	GOTLBO	RUN	JADE	PGJAYA	EMRUN	CGSCA
R	0.786	0.451	0.612	0.717	0.918	0.576	0.991	0.517
RMSE	2.216	3.363	2.855	2.543	1.386	2.915	0.475	3.328
MAPE	6385.878	11,701.750	12,167.122	9225.505	5000.666	12,235.108	415.259	10,848.074
NSE	0.600	0.078	0.336	0.473	0.843	0.307	0.982	0.097
IA	0.885	0.682	0.771	0.846	0.957	0.746	0.995	0.731
MdAE	0.760	1.671	1.502	0.984	0.397	1.366	0.021	1.687
MaxAE	8.126	17.090	11.393	10.515	4.191	11.308	1.764	15.510
WHD	39.234	44.415	45.098	43.174	30.731	44.827	10.962	43.076

lists the results of all optimization algorithms on the basis of the eight metrics. Among the methods tested, the EMRUN method, with the largest amount for R (0.991), $I_A$ (0.995), and $NSE$ (0.982) and the smallest value for $RMSE$ (0.475), $MAPE$ (415.26), $MdAE$ (0.021), $MaxAE$ (1.764), and $WHD$ (10.962), outperforms the ADEPSO, RUN, ATLDE, GOTLBO, JADE, PGJAYA, and CGSCA in deriving the optimal release from the six-reservoir problem. Consequently, these results demonstrate that the proposed method, with the RK-SM, IESQ, CMA-ES, and SQP operators, has a superior performance in optimally solving the six-reservoir problem.

Specifically, the Taylor diagram [57], a constructive approach, is used in this investigation to demonstrate how well the method performs compared with the other methods. This diagram shows that the closest point to the “reference” corresponds to the technique with the most outstanding performance. Taylor diagrams displayed in Figure 8 again demonstrate that the EMRUN algorithm has the most remarkable overall performance, followed by the JADE, ADEPSO, RUN, GOTLBO, PGJAYA, CGSCA, and ATLDE algorithms, in order.

5.2. Power-Generated Analysis

Figure 9 displays the power generated by all the optimization methods for six reservoirs. This figure is used to indicate the ability of the proposed EMRUN algorithm to maximize the power from the six-reservoir problem. According to the figure, the EMRUN algorithm can generate more power from all reservoirs except reservoir 3. Figure 10 depicts the monthly produced power that has been optimized using all optimization techniques. It can be seen in the figure that the quantity of power obtained by the EMRUN results in a remarkable agreement with the global solution. Consequently, these findings again demonstrate that the proposed method has high efficiency in managing the hydropower multireservoir system.

6. Conclusions

The optimization of water engineering systems is an essential study field for engineers. Water engineering problems are some of the most challenging engineering problems owing to their nonlinear, nonconvex, and high-dimensional natures. Hence, developing and creating an effective optimization method to solve these problems is an inevitable task. Therefore, this study effectively developed EMRUN to tackle water engineering problems. The innovation of the EMRUN lies mainly in improving the exploration phase by using the RK-SM and CMA-ES methods and in boosting the exploitation phase by utilizing the IESQ and SQP methods. In addition, adaptive parameters (i.e., $A{C}^{\prime }$, $\psi$, $p{c}_n$, and $g_n$) were implemented to improve the stability of these two phases. The advantage of EMRUN was initially shown using CEC-17 benchmark functions with various properties. The comparisons with the ADEPSO, JADE, GOTLBO, PGJAYA, RUN, CGSCA, and ATLDE verified EMRUN’s superior efficiency in exploration, exploitation, and local optimal avoidance.

Furthermore, an eight-parameter Muskingum model, called the generalized nonlinear Muskingum (GNLM4), was used to assess the proposed EMRUN algorithm to predict flood routing. In this work, two separate case studies were employed to assess the performance of the suggested algorithm. The two case studies were the Wilson flood (WF) and the Viessman and Lewis flood (VLF). For the WF, the findings indicated that EMRUN beat all other methods and was able to decrease the SSD by between 6% and 85%. In addition, in the case of VLF, the proposed EMRUN can reduce the SSD from 12% to 38% compared with the other methods. EMRUN findings also successfully predicted peak discharge, with a quantity remarkably similar to what was observed.

Finally, the proposed EMRUN is appraised to optimize the operation of a hydropower six-reservoir system. The fitness function amounts for EMRUN, the ADEPSO, JADE, ATLDE, RUN, GOTLBO, PGJAYA, and the CGSCA for this problem were 25.50, 24.62, 24.96, 22.12, 24.08, 23.08, 2280, and 22.37, respectively. According to the results, the EMRUN algorithm can more precisely and reliably optimize the six-reservoir system than the other algorithms can. Overall, the EMRUN method is superior to the other optimization methods because it relies on its effective operators (RK-SM, IESQ, CMA-ES, and SQP) and adaptive parameters, making EMRUN an effective, comparatively user-friendly optimization method with a high potential to optimize various types of water engineering problems. In practice, EMRUN seems to be a viable approach. In the future, we will use advanced-learning algorithms such as oppositional learning and orthogonal learning to boost EMRUN performance even more and then use them to solve the other water-resource challenges.