Comparing the Efficiency of Particle Swarm and Harmony Search Algorithms in Optimizing the Muskingum-Cunge Model

1. Introduction:

Flood routing in channels and rivers is of paramount importance and is utilized in nearly all river engineering studies for flood control, involving unsteady flow computations[1-3]. Typically, during flood events, the main river channel fills, and floodwaters inundate floodplains. Under these conditions, the hydraulics of flow become highly complex due to the abrupt changes in flow velocity within the main channel and floodplains. The continuity and momentum equations, which form the basis of unsteady flow models, have limitations in such scenarios. Even using equations like the Saint-Venant equation, a cornerstone of flood routing, can lead to inaccuracies. These discrepancies may arise from model parameters that are sensitive to the basin’s geography and flow conditions[4-6]. Various mathematical methods have been proposed to refine these equations. Most of these methods introduce additional variables and parameters, increasing model complexity and computational time while improving accuracy. Simplified methods with minimal data requirements can be more practical. Unlike hydraulic models that employ numerical methods for flood routing and involve complex unsteady flow equations, hydrological models utilize a spatially lumped continuity equation and a storage equation for flood routing, requiring less data for flood prediction[7,8]. One such method is the Muskingum-Cunge model, renowned for its simplicity and accuracy, making it one of the most widely used hydrological routing models. However, determining optimal values for its parameters, especially the weighting coefficient, remains challenging.

1.1. Flood Routing Models and the Position of the Muskingum-Cunge Model

Flood routing is a critical aspect of hydrological studies and water resources management. Various models are employed to simulate and predict the behavior of flood flows in rivers and channels. These models can be broadly categorized into hydraulic and hydrological models. Hydraulic models are based on the numerical solution of unsteady flow differential equations, such as the Saint-Venant equations or continuity and momentum equations. Examples include the dynamic wave model, kinematic wave model, and diffusion wave model. While offering high computational accuracy, these models are complex, require extensive input data, and have long computation times [9,10]. In contrast, hydrological flood routing models utilize simpler and more approximate forms of flow equations. Notable examples include the Muskingum model, Clark model, and Muskingum-Cunge model. Hydrological models are characterized by their computational efficiency, limited data requirements, and simplicity. However, their accuracy is generally lower than that of hydraulic models [11,12]. Among hydrological flood routing models, the Muskingum-Cunge model holds a special place. This model, an extension of the original Muskingum model, combines the continuity equation with a simple storage equation to simulate flood flow over a time period. The Muskingum-Cunge model requires only two primary parameters, namely the time lag and weighting coefficient, which can be determined through calibration with observed data. Its simplicity and relatively good accuracy have made it one of the most widely used tools for hydrological routing.

1.2. Calibration Methods for Hydrological Models

Calibration of hydrological models, including the Muskingum-Cunge model, is a critical step in their application. Calibration involves determining optimal parameter values by comparing model outputs with observed data. For the Muskingum-Cunge model, the time lag and weighting coefficient must be adjusted to achieve the best fit between the simulated and observed hydrographs. Calibration methods can be broadly classified into manual and automated approaches. Manual calibration involves manually adjusting parameter values and comparing the output with real-world data. This is a time-consuming and tedious process, especially for models with numerous parameters. Automated methods, on the other hand, seek to find a set of parameters that minimizes an objective function, such as the sum of squared errors, thereby achieving the best fit with observed data.Various automated calibration methods exist, each with its own advantages and disadvantages. One popular approach involves the use of evolutionary optimization algorithms, which have gained significant attention in recent years [12].

1.3. Application of Optimization Algorithms in Solving Hydrological Problems

Evolutionary optimization algorithms, such as genetic algorithms, particle swarm optimization, harmony search, and ant colony optimization, have demonstrated remarkable capabilities in solving complex multi-variable optimization problems. Consequently, their application in hydrological problems, including model calibration, parameter estimation, and water resources systems optimization, has garnered significant attention from researchers. Flood routing is a critical issue in the design and evaluation of hydraulic structures, water resources planning, and forecasting river flow and maximum water levels. Since most rivers lack continuous monitoring stations, methods are needed to estimate flow and water levels at specific locations. Flood routing is therefore used to assess flood wave profiles along a river and gain a better understanding of flood flow characteristics. Moreover, flood management can be viewed as risk management, involving measures to reduce the risk of floods or mitigate their impacts, such as minimizing damage to infrastructure, residential areas, and loss of life [13]. Generally, flood routing methods can be categorized into hydraulic and hydrological methods[14]. Hydraulic methods are based on the solution of differential continuity and momentum equations. While offering high accuracy, hydraulic methods require data such as flow depth, velocity, roughness coefficient, cross-sectional area, and bed slope [15]. Collecting such data is time-consuming and costly. However, hydraulic methods are relatively simple and reasonably accurate, requiring only one or more simultaneous input or output hydrographs from a gauging station [13]. The Muskingum method is one of the most well-known and widely used hydrological flood routing techniques. In this method, model parameters (decision variables) are estimated through optimization methods based on historical input-output hydrograph data. Then, the Muskingum routing equations are executed to determine an output hydrograph corresponding to an input hydrograph [16]. The Muskingum method can be applied to river routing in both linear and nonlinear forms. In the linear model, a linear relationship between storage and discharge is assumed, and subsequent calculations are based on this assumption. The linear Muskingum model has been used by various researchers [17,18]. Since the storage-discharge relationship is generally nonlinear[19], the linear Muskingum model may not provide sufficient accuracy in simulating floods. Therefore, Gill (1978) introduced a nonlinear version of the Muskingum method. The existence of a nonlinear relationship between storage and discharge in a river system is closer to reality and seems more logical. Consequently, the nonlinear Muskingum method has gained more attention from researchers in recent years [14,20,21]. Comparisons between the Muskingum-Cunge and kinematic wave methods have shown that, although the kinematic wave method utilizes more accurate data, the Muskingum-Cunge method is more efficient for flood routing. This efficiency is attributed to its lower data requirements and the generation of hydrographs that closely resemble real-world data compared to the kinematic wave method[22]. A four-parameter nonlinear Muskingum model has also been proposed, demonstrating more accurate predictions of flood discharge compared to the three-parameter model [23]. Moradi et al. (2008) applied both the Muskingum and Muskingum-Cunge methods to a section of the Liqvan River and compared the results for floods with return periods of 2, 10, 20, 50, and 100 years. The results revealed differences between the routing results of the two models and actual values. Fotuhi and Maghrebi (2010) compared dynamic and Muskingum-Cunge routing in a channel study and found minimal differences between the two methods. However, they cautioned against using the Muskingum-Cunge method in low-gradient channels. Norouzi and Bazargan (2020) investigated the accuracy of the Muskingum-Cunge method between the Malthan and Ahvaz gauging stations along the Karun River. Their findings indicated that the parameters of the Muskingum-Cunge method performed better when using three different values instead of a single constant value. This study aims to develop a model for the Balochistan region, prone to seasonal river flooding during monsoon seasons, using well-established algorithms and an optimization process. Given the limited data availability from flow gauging stations in Balochistan’s rivers, the Muskingum-Cunge model, which requires minimal data, is employed. Furthermore, to address the challenge of parameter estimation in the routing process, metaheuristic algorithms such as Particle Swarm Optimization (PSO) and Harmony Search Algorithm (HSA) will be utilized in conjunction with a three-parameter Muskingum model to optimize the routing in eight flood-prone rivers within the region.

2. Materials and Methods

2.1. Study Area

Sistan and Baluchestan province, spanning 180,726 square kilometers, is located in southeastern Iran, between 25°3’ and 31°27’ north latitude and 58°50’ and 63°21’ east longitude. The province comprises two main regions: northern Sistan and southern Baluchestan. It is divided into 15 counties (Zahedan, Mirjaveh, Khash, Saravan, Iranshahr, Nik Shahr, Chabahar, Sarbaz, Konarak, Delgan, Mehrestan, Sib and Soran, Bampoor, Qasr-e Qand, and Fanuj) (Statistical Center of Iran, 2020). This research focuses on the Baluchestan region, utilizing data collected from hydrometric stations along the Rapch, Sarbaz, Kehir, Siahjangal, Bahu, Kajou, Nikshahr, and Kamb rivers to conduct flood routing using the Muskingum-Cunge method (Figure 1). The characteristics of the studied rivers and their basins are presented in Table 1. These data will be used to construct inflow hydrographs for the river basins in the study area. Given the region’s vast geographical extent, it exhibits a diverse climate. Climate zoning indicates that the counties of Nikshahr, Qasr-e Qand, Mirjaveh, Konarak, Chabahar, Sarbaz, and Iranshahr have arid hot desert, semi-arid temperate, arid hot coastal, and extremely arid hot desert climates, respectively (Moradi et al., 2020). The region is characterized by high average temperatures (22.6°C) and low temperature fluctuations. Due to low annual precipitation (139.8 mm) and the absence of mountainous snow sources, most river flows are seasonal and temporary. In many parts of Baluchestan, limited groundwater resources serve as the primary water supply [12].

2.2. The Muskingum-Cunge Method

The Muskingum method, initially proposed in 1938 based on studies of the Muskingum River in the United States, is a linear model founded on the continuity equation (McCarthy, 1939). The model is represented by the following equations:

$${\displaystyle \frac{d{s}_t}{dt}}=I_t-O_t$$

(1)

$$S_t=K\left[x{I}_t+(1+x)O_t\right]$$

(2)

where:

S_t is the storage volume of the river (m³), I_t is the inflow rate (m³/s), Q_t is the outflow rate (m³/s)

K is the storage time constant (s), varying between 0 and 30 x is the weighting factor, ranging from 0 to 0.5

In 1969, Cunge demonstrated that this equation is analogous to the advection-dispersion equation and that its results closely resemble those obtained from the kinematic wave method. By modifying the kinematic wave equation and aligning physical diffusion with numerical diffusion, Cunge refined the Muskingum equation. The resulting modified Muskingum-Cunge equation, capable of considering lateral inflows, can be calculated based on the river’s physical characteristics (Cunge, 1969). This method has been widely used for flood routing in various regions. Assuming that inertia and acceleration can be neglected compared to the bottom slope in natural channels, Cunge proved that the advection-dispersion equation is similar to the Muskingum equation. Based on this, the Muskingum-Cunge equation for flood routing can be defined as follows (Cunge, 1969):

$$Q_{i+1}^{j+1}={\alpha }_1{Q}_i^{j+1}+{\alpha }_2{Q}_i^j+{\alpha }_3{Q}_{i+1}^j$$

(3)

where the coefficients ${\alpha }_1,{\alpha }_2,$ and ${\alpha }_3$can be calculated using equations (4) to (6) (Chow, 1959):

$${\alpha }_1=\frac{\Delta t-2KX}{2K(1-X)+\Delta t}$$

(4)

$${\alpha }_2=\frac{\Delta t+2KX}{2K(1-X)+\Delta t}$$

(5)

$${\alpha }_3=\frac{2K(1-X)-\Delta t}{2K(1-X)+\Delta t}$$

(6)

In these equations, Δt is the time step, Δx is the spatial step, and K and x are routing parameters. It should be noted that the sum of ${\alpha }_1,{\alpha }_2,$ and ${\alpha }_3$is equal to 1, and both K and Δx have units of time. In the above equations, K is the storage constant, and x is a factor that represents the relative influence of inflow on storage and is dimensionless. Cunge showed that when K and x are assumed to be constant, equation (3) is an approximate solution to the kinematic wave equations. Additionally, this equation is also an approximate solution to the diffusion equation in a channel. Therefore, the values of the routing parameters K and x can be calculated using the following equations:

$$K={\displaystyle \frac{\Delta X}{C_K}}={\displaystyle \frac{\Delta X}{\frac{dQ}{dA}}}$$

(7)

$$X={\displaystyle \frac{1}{2}}(1-{\displaystyle \frac{Q}{B{C}_{K}S\Delta x}})$$

(8)

where $C_K$ is the wave celerity related to the flood discharge, B is the water surface width in the river or channel, Q is the reference discharge, and Δx is the spatial step. The reference discharge is the discharge used to accurately estimate the routing parameters (Akbari et al., 2010). Equation (8) shows that the value of x is directly related to the reference discharge, which is a function of the inflow. Since the inflow varies during a flood event, the reference discharge also varies. Therefore, in the constant parameter Muskingum-Cunge method, the reference discharge is considered constant in all computational cells. In this study, the following equation is used to estimate the reference discharge (Akan et al., 2006; Wilson & Ruffini, 1988):

$$Q=Q_b+0.5\left(Q_{Pi}-Q_b\right)$$

(9)

where Q_b and Q_pi represent the minimum and maximum inflow rates, respectively.

It is observed that estimating the routing parameters K and x is the most critical part of flood routing estimation, and traditional manual methods are time-consuming. Therefore, modern methods based on artificial intelligence should be used to estimate these parameters, as they are more accurate and time-efficient. In this study, two optimization algorithms, particle swarm optimization (PSO) and harmony search (HS), will be used to estimate these parameters and find the best solution. The programming of these machine learning algorithms will be done in the MATLAB programming environment.

2.2. Particle Swarm Optimization (PSO)

PSO is an evolutionary algorithm introduced in 1995. Based on population-based search, similar to genetic algorithms, ant colony optimization, and bee colony optimization, PSO is a nature-inspired algorithm that mimics the collective intelligence and social behavior of birds and fish [24,25]. Due to its simple structure, few controllable parameters, rapid convergence, and high computational efficiency, PSO is widely used in optimization problems. Although primarily developed for optimization problems with continuous variables, PSO has been successfully adapted for discrete variable problems. In PSO, each solution is considered a bird in a flock and is termed a particle. These particles exhibit both individual intelligence and social behavior, coordinating their movement towards a target (Montalvo et al., 2008). The process begins with a swarm of particles, each representing a randomly generated solution to the problem. The swarm then iteratively searches for the optimal solution. Each particle i is defined by three vectors: its current position Xi, its best previous position Yi, and its velocity Vi.

$$X_i=(x_{i1},x_{i2}\mathrm{.....}{x}_{in})\mathrm{current}\mathrm{position}$$

(10)

$$|Y_i=(y_{i1},y_{i2}\mathrm{.....}{y}_{in})\mathrm{Best}\mathrm{previous}\mathrm{position}$$

(11)

$$V_i=(v_{i1},v_{i2}\mathrm{.....}{v}_{in})\mathrm{Flight}\mathrm{velocity}$$

(12)

It is assumed that the simulated bird flock communicates while flying, and each bird looks towards a specific direction, which is the best position achieved (Yi) in the current state. Subsequently, when the birds communicate, the bird in the best position is identified. Through coordination, each bird moves towards the best bird with a velocity that depends on its current velocity. Therefore, each bird explores the search space from its current local position, and this process repeats until the bird potentially reaches the desired position. Birds learn through their own experience (local search) and the experience of their peers (global search). In each iteration, the particle with the best momentary solution to the problem is identified. The position of this particle is then used in calculating the new position for each particle in the swarm. This calculation is based on the following equation:

$$X_{i,d}^{n+1}=X_{i,d}^n+V_{i,d}^{n+1}$$

(13)

The best position achieved by each particle in each swarm is denoted by p-best, and the best position achieved in the entire population is denoted by g-best. Ultimately, the coordinated population moving in one direction is represented by the following equation:

$$V_{i,d}^{n+1}=\omega V_{i,d}^n+c_1{r}_1^n(pbes{t}_{i,d}^n-X_{i,d}^n)+c_2{r}_2^n(gbes{t}_{i,d}^n-X_{i,d}^n)$$

(14)

where d=1, 2, …, D, i=1, 2, …, N, N is the population size, ω is the inertia weight, c₁ and c₂ are positive constant coefficients that refer to cognitive and social parameters (learning factors), respectively, r1 and r2 are random numbers between 0 and 1, and n=1, 2, … denotes the number of iterations [26].

Furthermore, inertia weight is an index to ensure the convergence behavior of PSO, which is used to control the effect of the history of velocities on the current velocity. This coefficient may change from one iteration to the next. Since this coefficient controls the balance between global and local search, it is suggested to decrease linearly with time, usually in a way that global search is prioritized at the beginning and local search is prioritized with each subsequent iteration[27]. To determine the inertia weight, the following equation is used:

$$w=w_{\max }-{\displaystyle \frac{(w_{\max }-w_{\min })\times n}{iTe{r}_{\max }}}$$

(15)

where ωmax is the initial inertia, ωmin is the final inertia, iTer_max is the maximum number of iterations, and n is the current iteration number. According to experimental studies, the best values for ωmax and ωmin are 0.9 and 0.4, respectively [28]. Although the constant coefficients c1 and c2 do not have a decisive role in the convergence of PSO, by determining appropriate initial values, it is possible to increase the convergence speed and improve the local solution. The suggested initial values for these coefficients also show that choosing a larger cognitive parameter (c1) compared to the social parameter (c2), with the condition c1 + c2 ≤ 4, leads to better results. The flowchart of the particle swarm optimization algorithm is shown in Figure 2

3.2. Harmony Search Algorithm: A Musical Inspiration

The Harmony Search (HS) metaheuristic algorithm draws inspiration from the musical process of searching for a state of “perfect” harmony, akin to jazz improvisation [29]. Jazz improvisation seeks the optimal state (a superb harmony) as determined by aesthetic evaluation, just as an optimization algorithm seeks the optimal state (global optimum) as determined by an objective function evaluation. Aesthetic evaluation is performed by the set of notes played by each instrument, similar to how the objective function is evaluated by a set of values determined by each decision variable. The quality of the harmony improves with consistent practice, just as the quality of the solution improves with iterations.

Key Advantages of HS:

Ability to handle discrete variables: Unlike many numerical optimization methods, HS can effectively handle discrete variables, reducing the likelihood of getting trapped in local optima [30,31].

Steps Involved in the HS Algorithm:

Initialization: Set initial values for algorithm parameters such as Harmony Memory Size (HMS), Harmony Memory Considering Rate (HMCR), Pitch Adjusting Rate (PAR), and generate an initial harmony.

Initialization of Harmony Memory (HMS): Populate the HMS with a set of randomly generated harmony vectors.

Creating a new harmony using HMS: For each variable of a new harmony vector, three operators are employed: HMCR, PAR, and random selection.

-

HMCR: If a random number is less than HMCR, a value is selected from the HMS.

-

PAR: If a random number is less than PAR, the selected value is slightly adjusted.

-

Random selection: If neither HMCR nor PAR conditions are met, a random value within the variable’s bounds is selected.

Updating the Harmony Memory: Evaluate the new harmony using the objective function. If it’s better than the worst harmony in the HMS, replace the worst harmony with the new one.

Termination: Check if the termination criterion (e.g., maximum number of iterations) is met. If not, go back to step 3.

Here’s the translation and formulas for the evaluation metrics of machine learning algorithms:

3. Evaluation of Machine Learning Algorithm Performance

To evaluate the algorithms used, the following indices were employed: Coefficient of Residual Mass (CRM), Normalized Root Mean Square Error (NRMSE), Model Efficiency (EF), and Index of Agreement (d) [32,33].

1. Coefficient of Residual Mass (CRM):

CRM = (Σ Oi - Σ Pi) / Σ Oi

(16)

2. Normalized Root Mean Square Error (NRMSE):

NRMSE = √[(Σ (Pi - Oi)²) / n] / Ō

(17)

3. Model Efficiency (EF):

EF = 1 - [Σ (Oi - Pi)² / Σ (Oi - Ō)²]

(18)

4. Index of Agreement (d):

d = 1 - [Σ (Pi - Oi)² / Σ (|Pi - Ō| + |Oi - Ō|)²]

(19)

where:

Oi: Observed values, Pi: Predicted values, Ō: Mean of observed values, n: Number of observations

These metrics provide a comprehensive evaluation of the model’s performance, including its bias (CRM), accuracy (NRMSE), efficiency (EF), and overall agreement with observed data (d). NRMSE (Normalized Root Mean Square Error) represents the difference between simulated and observed data relative to the mean of observed data. For this index:

-

Values less than 10%: Excellent

-

Values between 10% and 20%: Good

-

Values between 20% and 30%: Fair

-

Values greater than 30%: Poor

The EF (Model Efficiency) index indicates the degree of alignment between observed and predicted plots with the 1:1 line. Its values range from 1 to negative infinity. An EF value of 1 corresponds to a perfect fit [34]. The d (Index of Agreement) ranges from 0 to 1:

-

1 indicates perfect agreement between simulated and observed data

-

0 indicates no agreement

-

For CRM (Coefficient of Residual Mass):

-

Positive values indicate underestimation of crop performance

-

Negative values indicate overestimation of crop performance

These metrics collectively provide a comprehensive evaluation of model performance, assessing accuracy, bias, and agreement between predicted and observed values. Finally, in addition to examining and evaluating the performance of the Muskingum-Cunge model for estimating the flood input hydrograph in the Baluchistan region, the routing parameters in this model are estimated using machine learning algorithms. By comparing these machine learning methods, the best approach for this region is determined for estimating the flood input hydrograph and its routing.

4. Results and Discussion

This study aims to find the optimal parameters (K and X) for the Muskingum-Cunge model using PSO and HS algorithms. The performance of these two algorithms in terms of convergence rate, accuracy, and objective function values is compared. The results of both algorithms, including convergence trends and optimal parameter values, are presented and discussed.

This Table 2. presents initial estimates for parameters X and K for eight rivers, derived using Particle Swarm Optimization (PSO) and Harmony Search (HS) algorithms. The results show similar values between PSO and HS for most rivers, suggesting consistency across methods. Parameter X values are relatively stable across all rivers, ranging from about 0.16 to 0.43. In contrast, parameter K shows greater variation, with values from 2.85 to 29.11.

Figure 4, presents a Comparison of PSO and HS Optimization Models with Actual Flood Discharge Data This analysis compares the results from Particle Swarm Optimization (PSO) and Harmony Search (HS) optimization models with actual flood discharge data across multiple stations. The performance of these models was evaluated using correlation coefficients (R²) and visual comparison of simulated hydrographs with actual flood discharge data.

Kajoo Station: Both PSO and HS models simulated the flood discharge with high accuracy. Peak discharge values were 231.60, 234.07, and 236.52 m³/s for actual, PSO, and HS models respectively. Correlation coefficients were exceptionally high at 0.9458 (PSO) and 0.9451 (HS), indicating very strong correlation between simulated and actual data.

Bahu Station: Both algorithms successfully simulated the flood output hydrograph with acceptable accuracy. The simulated graphs closely matched the actual flood discharge graph. Correlation coefficients were 0.6879 for PSO and 0.6558 for HS, which are close to 0.7, indicating acceptable correlation between simulated and actual data.

Nikshahr Station: The algorithms simulated the flood hydrograph with acceptable accuracy. Correlation coefficients were 0.8031 for PSO and 0.7790 for HS, showing moderate to good correlation.

Rapch Station: Simulated hydrographs closely matched actual flood discharge. Correlation coefficients of 0.7834 (PSO) and 0.7747 (HS) demonstrated moderate to good correlation.

Sarbaz Station: Both algorithms provided acceptable accuracy, with correlation coefficients of 0.7651 (PSO) and 0.755 (HS), showing moderate to good correlation.

Siah Jangal Station: This station showed the highest correlation among all stations, with coefficients of 0.9645 (PSO) and 0.9631 (HS), indicating very high agreement between simulated and actual data.

Kahir Station: Both algorithms achieved highly accurate simulation, with correlation coefficients of 0.8568 (PSO) and 0.8495 (HS), showing very good correlation.

Kamb Station: While overall simulation was acceptable, both algorithms struggled to accurately simulate flood discharge during certain time periods, resulting in significant differences between simulated and actual data in these instances.

In conclusion, both PSO and HS optimization algorithms demonstrated good overall performance in simulating flood hydrographs across the studied stations. The Siah Jangal station showed the highest correlation, while the Bahu and Kamb stations presented some challenges in accurately simulating certain aspects of the flood discharge. These results suggest that both PSO and HS are effective tools for flood discharge simulation, although their performance may vary depending on the specific characteristics of each station and time period.

Figure 5 illustrates the charts of inlet flow and outlet flow with PSO and HS models, optimizing the Muskingum-Cunge coefficients for various stations. Kajo Station: The chart shows the input and output discharge using coefficients optimized by PSO and HS models. Peak outflow discharge values for PSO and HS models were 172.45 and 174.10 cubic meters per second, respectively, indicating acceptable accuracy for both models in simulating the outflow discharge. Bahu Station: Both PSO and HS algorithms, by optimizing Muskingum-Cunge coefficients, have simulated the flood outflow discharge with acceptable accuracy. The simulated outflow discharge graphs for both algorithms are close to each other. NikShahr Station: The optimization of Muskingum-Cunge coefficients by both PSO and HS algorithms has resulted in acceptable accuracy in simulating flood outflow discharge. The simulated outflow discharge graphs for both algorithms are similar. Rapch Station: Both PSO and HS algorithms have successfully optimized Muskingum-Cunge coefficients, leading to acceptable accuracy in simulating flood outflow discharge. The simulated outflow discharge graphs for both algorithms are close to each other. Sarbaz Station: The optimization of Muskingum-Cunge coefficients by PSO and HS algorithms has produced acceptable accuracy in simulating flood outflow discharge. The simulated outflow discharge graphs for both algorithms show close agreement. Siyah Jangal Station: Both PSO and HS algorithms have optimized Muskingum-Cunge coefficients very effectively, resulting in highly accurate simulation of flood outflow discharge. The simulated outflow discharge graphs for the two algorithms show very close alignment with each other. Kahir Station: The optimization of Muskingum-Cunge coefficients by PSO and HS algorithms has led to good accuracy in simulating flood outflow discharge. The simulated outflow discharge graphs for the two algorithms show good agreement, indicating similar performance. Kamb Station: Both PSO and HS algorithms have successfully optimized Muskingum-Cunge coefficients, resulting in good accuracy in simulating flood outflow discharge. The simulated outflow discharge graphs for the two algorithms show good agreement, demonstrating good and similar performance. Overall, these results confirm that both PSO and HS algorithms perform well in optimizing Muskingum-Cunge coefficients across various stations, producing accurate simulations of flood outflow discharge.

Heat maps in Figure 6 compare the performance of PSO and HS algorithms in optimizing the Muskingum-Cunge model coefficients at different stations based on CRM, Ef, d, and NRMSE indices. Brighter colors indicate better performance. Generally, PSO outperforms HS in terms of CRM, Ef, and NRMSE at most stations, suggesting its superiority in simulating flood hydrographs. However, there’s little difference between the two algorithms in terms of the ‘d’ index. At stations like Kajo, NikShahr, Rapch, Sarbaz, Siyah Jangal, and Kahir, PSO consistently provides more accurate simulations. At Kamb, PSO significantly outperforms HS across all indices, particularly in terms of Ef (PSO: -2.69 vs HS: -14.10) and NRMSE (PSO: 0.24 vs HS: 0.49). The performance difference is most pronounced at this station, with PSO showing markedly better results in simulating flood hydrographs.

Table 3. presents statistical analysis results for four indices (CRM, Ef, d, NRMSE) comparing PSO and HS models at various stations. PSO consistently performed as well or better than HS across all locations. Notably, PSO significantly outperformed HS at Bahu and Kamb stations, with substantial improvements in Ef and CRM indices. While both models performed well at Siyah Jangal, PSO showed marginally better results. At Kajo and NikShahr, PSO demonstrated moderate performance advantages across the indices

5. Conclusion

This study compared the performance of Particle Swarm Optimization (PSO) and Harmony Search (HS) algorithms in optimizing Muskingum-Cunge model parameters for flood routing across multiple stations. Results consistently demonstrated PSO’s superiority over HS in simulating flood hydrographs. While both algorithms achieved acceptable accuracy at most stations, as evidenced by correlation coefficients ranging from 0.6558 to 0.9645, PSO exhibited consistently better performance across various evaluation metrics (CRM, Ef, NRMSE). This superiority was particularly pronounced at stations with complex hydrological conditions, such as Bahu and Kamb. Both PSO and HS effectively optimized Muskingum-Cunge model parameters, leading to accurate simulations of flood outflow discharge. However, PSO’s ability to explore the solution space more efficiently resulted in superior parameter estimates and consequently improved hydrograph predictions. Heat map analysis and statistical comparisons further reinforced PSO’s dominance. While both algorithms captured the general trends in flood hydrographs, PSO exhibited a more precise representation of peak flows and hydrograph shapes.While both PSO and HS are valuable tools for flood modeling, the findings of this study strongly advocate for PSO as the preferred optimization algorithm for the Muskingum-Cunge model. Future research should focus on understanding the underlying reasons for PSO’s superior performance and exploring hybrid approaches that combine the strengths of both algorithms.