Multi-objective optimization of stepped spillway and stilling basin dimensions

In the present paper, an optimization procedure is proposed for stepped spillway design dimensions, which leads to maximum energy dissipation rate and minimum construction cost, considering independently the chute cost and stilling basin cost. Three independent objective functions are thus simultaneously satisfied. The procedure involves four main tools: the multi-objective particle swarm optimization method (MOPSO) to find Pareto solutions in one run, the K-means clustering algorithm to reduce the size of the obtained non-dominated solutions, the pseudo-weight vector approach (PWV) to facilitate the decision making and to select some adequate solutions, and finally, CFD simulations to analyze the retained optimal solutions. The suitability of the proposed procedure is tested through an example of application. As results, a set of 20 solutions with different satisfaction levels is found and compared to existing solutions. A multi-objective optimization problem may have many different solutions; the originality of the present work lies in the proposed procedure, which explores several possible solutions and reduces their number to give help to the decision making. Furthermore, an approximate expression of spillway total cost is also derived as a function of flow energy dissipation rate.


INTRODUCTION
Construction of dams serves multiple purposes, including water supply, river flow control, flood control, power generation, etc. To avoid damage or collapse of dams by overflowing, an adequate ancillary structure is required for efficient flood release. The stepped spillway is one of the most important of these structures, which stands for its high energy dissipation capacity (Chamani & Rajaratnam 1994). The design of stepped spillways dates back to the early times. Different civilizations have contributed, for more than 3,000 years, to the development of the art of building dams and stepped spillways (Chanson 2001).
Over the centuries, stepped chutes were a common conception of hydraulic structures, but at the beginning of the 20th century, breakthroughs in the understanding of hydraulic jumps favored the design of stilling basins to the detriment of these structures (Chanson 2001). However, with the introduction of roller compacted concrete (RCC) and new construction technology, the design of stepped spillways has regained interest due to their advantages compared to other types of flood release structures (Chamani & Rajaratnam 1999). The main feature of stepped spillways is the high energy dissipation rate. The size of the stilling basin and therefore the whole construction cost are also reduced (Chanson 1994a). In addition, due to the large amount of air entrainment in the flow, the risk of cavitation is minimized and downstream bed erosion risk is reduced. It is well known that cavitation is an essential problem that should be avoided in designing all kinds of hydraulic structures and machinery (Al-Obaidi 2020a, 2020b, 2020cAl-Obaidi & Mishra 2020).

Stepped spillway design and modeling
In dam project, spillways are designed to release excess discharge with safety and sufficient capacity to avoid any overflowing. Their type and dimensions are determined with respect to cost consideration, hydrology conditions, dam height, topographic features, and foundation geology. Nevertheless, flow energy dissipation within a stepped spillways chute remains one of the most important parameter in the design procedure. Reducing flow energy allows reduction of cavity risk and erosion effects in the walls and the downstream river bed. In addition, for more dissipation and to ensure more protection of the downstream dam bottom, a stilling basin is commonly designed according to the residual flow energy at the chute end. Consequently, three dimensions are to be determined: the spillway width b, the steps height h s (or their number, knowing the dam height H D ), and the basin length L B . The spillway slope f, which is the remaining design parameter of the structure, is often equal to the downstream slope of the dam. Otherwise, it is imposed by the topography conditions of the construction site, especially when a lateral spillway is adopted. Hydrologic modeling and assessment of the drained watershed lead to the maximum discharge Q to be released. This gives the unitary flow rate q ¼ Q=b to be used for the critical water depth evaluation h 3 c ¼ q 2 =g (the gravity constant g is taken as equal to 9.81 m/s 2 ). The critical height h c is a key parameter to distinguish the flow regimes, and then the adequate equations to be applied. Indeed, three flow regimes (nappe, transition and skimmed) may occur in stepped spillways depending on the geometrical characteristics, the wall roughness, and the flow rate (Rassaei & Rahbar 2014).
The amount of energy dissipation depends on the flow regime: it is more important in the nappe and in the transition regimes. However, due to flow instabilities, most stepped spillways are designed for maximum discharge outside the transition regime. There are several empirical formulas to evaluate the flow energy dissipation in stepped spillways, the choice of adequate ones to apply must be done with caution. In the present work, expressions proposed by Chanson (1994a) are considered for the two cases of nappe and turbulent flow regimes, and those proposed by Boes and Hager (Boes & Hager 2003) for the turbulent regime case. These later make a distinction between uniform and non-uniform turbulent flows.
The nappe flow regime appears for low flow rates, when h c =h s , 0:0916(tanf) À1:276 , with h s being the steps height. It is characterized by a series of waterfalls from one step to another followed by the formation of partially or fully developed hydraulic jumps. Energy dissipation, in this case, is given by Equation (1) in terms of the rate of residual energy E r at the outlet to the maximum energy E m available at the inlet (Chanson 1994a In the case of the turbulent flow regime, water falls in cascades like a coherent stream; this regime is initiated by large discharges when h c =h s . 1:057 À 0:465 tanf. The rate of the residual energy is given in Equation (2) (Chanson 1994a).
Expressions given by Boes & Hager (2003) for this flow regime are displayed in Equations (3) and (4) for the uniform and the non-uniform cases respectively.
In the above equations, f b is the coefficient factor of the wall roughness, and D h,w is the hydraulic diameter of the uniform flow.
Even though stepped spillways dissipate a large amount of energy before the discharge is returned to the downstream river, the residual energy, as evaluated at the chute bottom, should be dissipated at the downstream structure in order to avoid scouring. In the present work, a rectangular stilling basin is considered for this aim, in which the torrential regime is transformed to a fluvial regime by hydraulic jump. Its length, L B , is thus taken as greater than the jump length. It depends on the Froude number F r of the incident flow and is proportional to the water depth h p at the chute bottom (Durand et al. 1999).
for F r 1 12 h p for 1 , F r 1:7 X 1 (F r )X 5 (F r )h p for F r . 1:7 8 < : (5) h p is the solution of the cubic equation 2gh 3 p À 2gE r h 2 p þ q 2 ¼ 0, and X 1 (F r ) and X 5 (F r ) are empirical adjustment polynomials: The rate of residual energy evaluated by applying the two above empirical models is compared to experimental results taken from Salmasi & Özger (2014). For this purpose, 77 different cases are used. The predicted and the observed values are plotted in Figure 1, which shows that the empirical equations of Boes & Hager (2003) can predict the energy dissipation with more accuracy. The slope of the linear fitting is closer to one and computed errors are much smaller. Then, to get more insight in the comparison and to validate at the same time, the numerical modeling achieved in the present work, nine (9) and five (5) refined CFD models using Ansys-Fluent software are implemented and executed for the turbulent and the nappe flow cases, respectively. The obtained simulated results are shown in Figure 2 together with empirical predictions. It is recalled that the Boes and Hager model does not apply for the nappe flow case. It can be seen from this figure that the present numerical results and the empirical predictions are very close together. Hence, for the optimization design computations, the Chanson model is adopted for the nappe flow regime, while the Boes and Hager model is adopted for the turbulent flow regime. This model best matches experimental data and simulated numerical results.
These results are obtained for a spillway of 5.5 m height and 22:5 of mean slope with 6 steps, releasing a discharge Q ¼ 80 m 3 =s. It is also to be noticed that in the present work, more attention is paid to the turbulent flow case. The turbulence is accounted by using the standard k-1 model. The water free surface detection and tracking is achieved using air-water twophase flow mixture model. This model solves the conservation equations for the mixture, the volume fraction equations for the secondary phases, as well as algebraic expressions for the relative velocities.  (Chanson 1994a;Boes & Hager 2003) and measured ones are taken from Salmasi & Özger (2014).

Spillway optimization
In hydraulic computations, head loss is commonly used to quantify all energy dissipation occurring in a flow. However, in practical optimization applications, the rate of the flow residual energy E r =E m is preferably used as the objective function instead of the energy dissipation rate. This gives a minimization problem which is commonly handled by optimization methods. Equations (1), (3) and (4) are thus well written as a first objective function of the present multi-objective optimization problem. The two other functions are the costs of the spillway chute and the stilling basin. The construction cost of the spillway chute, denoted by C S , is assumed to be proportional to the concrete volume of all steps, and the construction costs of the stilling basin (C B ) depends on its length L B , on its wall thickness e, and on the width b. These quantities are expressed by: Figure 3(a) and 3(b) take into account the fact that the residual energy rate and the stilling basin construction cost are not defined for the case of transition flow regime. This regime is represented by the blank zone. The right zone corresponds to the nappe flow regime and the left, which is larger, represents the turbulent flow regime. Most of the optimal design dimensions belong to the latter. However, the chute construction cost C S (Equation (6)) is directly proportional to the product b h s , its isovalues are thus simply a set of hyperbolas well defined for all practical values of design dimensions (Figure 3(c)). It can be seen from Figure 3(a)-3(c), that the stepped spillway is as efficient in dissipating energy the greater its width, but the more expensive the stilling basin construction cost. The effect of the steps' height h s on both energy dissipation rate and  Corrected Proof basin cost is practically negligible for small values of b. This effect becomes perceptible only for great values of b, where it induced great increase of the cost for small gain of energy dissipation. On the contrary, both h s and b have the same effect on the spillway chute cost, and it is as significant as these quantities are small. These conflicting tendencies make it hard to predict the optimal design dimensions giving a satisfactory energy dissipation rate with acceptable construction cost.
Finding values of h s and b that minimize the three objective functions can be formulated as a nonlinear multi-objective optimization problem as follows: b max is the maximum width that the spillway can have. It is theoretically limited by the total dam length when the spillway is integrated to the dam. Otherwise, it depends on the site topography. h max s ¼ 1=2H D ; it means that at least there must be a minimum of two steps to consider the spillway as a stepped one. In addition to the two above constraints, to avoid the transition flow regime, only solutions satisfying (1:057 À 0:465 tanf) À1 h c h s 12:34 (tanf) 1:276 h c are accepted. The proposed procedure to solve this optimization problem follows four main steps, as summarized in Table 1 and illustrated by the flowchart in Figure 4. The details on the achieved tasks are given in the next sections.

Multi-objective particle swarm optimization
Particle swarm optimization is a heuristic search technique introduced in 1995 by Kennedy and Eberhart to solve optimization problem (Coello Coello & Lechuga 2002). This algorithm is driven by the social behavior of insect swarms and their coordinated movement. PSO is a population-based research algorithm in which particles change their position with time according to their own and their neighbor's experience.

Corrected Proof
Considering a d-dimension research space with n particles, each i th particle of the population (swarm) is described by its position X i ¼ (X i1 , X i2 , . . . , X id ) and its velocity V i ¼ (V i1 , V i2 , . . . , V id ), with i ¼ 1 Á Á Á n. In the present work d ¼ 2, and the number of optimal variables is 2 (h s and b). The best recorded performance of each particle and the global best performance of each d-dimension among the swarm are denoted by P best i ¼ (P best i1 , P best i2 , Á Á Á P best id ) and G best ¼ (G best 1 , G best 2 , Á Á Á G best d ), respectively.
The velocity and the position of each particle are updated using its previous position and its actual velocity. The new velocity for each particle takes into account its previous velocity, its best encountered position, and its neighbor's positions (Fallah-Mehdipour et al. 2011).
where j ¼ 1, 2, . . . d, and i ¼ 1, 2, . . . n. W is the inertial weight, c 1 and c 2 are positive constant parameters called acceleration coefficients, r 1 and r 2 are random values between 0 and 1. X t ij and V t ij , are, respectively, the velocity and the current position of variable member j of particle i at iteration t.
The success of PSO for both continuous nonlinear and discrete binary single objective optimization and its high speed of convergence attracted researchers to extend the use of this technique to other areas. It is recognized that Moore and Chapman proposed in 1999 the first extension of PSO strategy for solving multi-objective problems. Then other reported studies extended PSO to various multi-objective problems (Reyes-Sierra & Coello Coello 2006;Fallah-Mehdipour et al. 2011;Hojjati et al. 2018).
Before applying PSO in multi-objective problems, the algorithm must be modified. The major modification concerns the selection of P best and G best . In MOPSO, instead of a single solution a set of solutions is determined, also called a Pareto Corrected Proof optimal set. When solving single objective optimization problems, the leaders that each particle uses to update its position are completely determined once a neighborhood topology is established. In MOO, each particle might have a set of different leaders from which only one can be selected for position updating. Such sets of leaders are usually stored in different places from the swarm that is called an external archive (this is a repository in which the non-dominated solutions found so far are stored). The solutions contained in the external archive are used as leaders when the positions of particles of the swarm have to be updated. The contents of the external archive are also usually reported as the final output of the algorithm (Coello Coello & Lechuga 2002;Hojjati et al. 2018). The steps of the MOPSO algorithm are: 1. Input parameters and specify the lower and upper boundaries of each variable. 2. Initialization of swarm and archive: -Initialize the current position of each particle of the swarm randomly within the specified boundaries, and set the velocity to zero. -Evaluate each particle in the population.
-Initialize the set of leaders with non-dominated particles and store it in the repository (REP).
-Define each particle's coordinates according to the values of its objective function by generating hypercubes of the search space explored so far, and locate the particles using these hypercubes as a coordinate system. -Initialize the memory of each particle with a single local best for each particle and store it in the repository, this memory serves as a guide to move through the search space. 3. For each particle in the swarm: -Determine the best global position (G best ) for each particle i from the repository REP by randomly selecting a particle form within a hypercube. -Update the speed and the new position, within the specified boundaries, of each particle using Equation (8). 4. Update the archive of non-dominated solutions: -Evaluate each particle in the population.
-Update the contents of the repository REP together with the geographical representation of the particles within the hypercubes. -Update the memory of each particle by using Pareto dominance. 5. Repeat steps 3 and 4 till the maximum number of iterations is reached.

Solution clustering and weighting
MOPSO, as the majority multi-objective optimization algorithms, is capable of finding multiple and diverse Pareto optimal solutions in one run, in which one objective cannot be improved without sacrificing other objectives. A large population size is thus necessary for a successful run and convergence of MOPSO. The population size cannot be chosen according to a small desired number of non-dominated solutions. In practical applications, a user is interested in a handful of solutions. For this reason, clustering algorithms are used to reduce the size of the obtained non-dominated set of solutions to a few representative ones. The most popular clustering method is the K-means algorithm (Santhanam & Velmurugan 2010). The steps of this algorithm are as follows (Deb 2001): 1. Initialize the number K of desired clusters and randomly generate K centroids: C 1 , C 2 , Á Á Á, C K . 2. Assign each particle, of position X i , from the total non-dominated N solutions, to a group that has the nearest centroid C k . This is achieved by measuring the Euclidean distances: 3. Update the centroids by taking the average position of each cluster. 4. Repeat steps 3 and 4 until the convergence criterion is met: centroids don't change with iteration.
After clustering, the number of solutions can be greatly reduced without loss of information. To facilitate decision making and to understand how each solution can influence the decision, the pseudo-weight vector approach may also be employed. The main idea of this approach is to calculate the relative distance of the solution from a maximum value in each objective function and for each obtained solution.
For minimization problem, the pseudo-weight is calculated for each objective function by the following equation (Deb 2001): Calculation of weights makes it easy for an end user to take a decision and to select a convenient or a preferred combination. In general, a solution with close weights values is a good compromise that satisfies all the considered objective functions.

APPLICATION AND RESULTS DISCUSSION
The suitability of the proposed optimization procedure is shown herein through an example of a stepped spillway of total height H D ¼ 5:5 m and mean slope f ¼ 22:5 . The design total discharge is Q ¼ 80 m 3 =s. This example was first introduced and solved using GA by Sharifi et al. (2005) and Haddad et al. (2005). It was then reproduced and solved using HBMO by Haddad et al. (2010). It was found that this well detailed and easily reproducible example constitutes a good comparison test for the proposed procedure.
The aim of the present application is to determine the best combinations of the spillway width b and the step height h s that minimize both the residual energy ratio and the total construction cost including, independently, the chute and the stilling basin costs. The thickness of the basin walls is first assumed as unity (e ¼ 1 m), and then complementary results are given for other values.
The procedure is implemented in the Matlab programming environment, as a set of scripts including objective functions definition, MOPSO algorithm, K-means clustering, pseudo-weights vector approach, and several plotting codes. A set of 500 Pareto optimal solutions was evaluated using 1000 MOPSO iterations. These numbers of particles and iterations were adopted after several numerical tests. These great values ensured the convergence of the solution. The following values W ¼ 0:5, c 1 ¼ 1, c 2 ¼ 2, r 2 ¼ r 1 ¼ 0:4 are used for the MOPSO parameters.
Solving this three-objective optimization problem led the non-dominated solutions to be obtained in one run. The great number of 500 Pareto solutions is convenient only for effectiveness of MOPSO searching, it is not adequate for decision making. It is thus reduced to 20 dominated solutions by clustering using K-means algorithm. Then they are assigned weights by applying a pseudo-weights vector approach to quantify the predominance of the objective functions over each other. 3Dspatial representation of the obtained Pareto and clustered solutions is given in Figure 5(a) and 5(b), respectively. It can be seen that the overall dense shape of the 500 solutions is globally well represented by a small number of particles that are well positioned in the (f 1 , f 2 , f 3 )-pseudo 3D space. Table 2 reports values of these clustered solutions for the three objective functions and their corresponding weights. The values found for the geometric dimensions (b, h s and L B ) are also given. This kind of table is of great help for the decision making regarding the solution to be retained for real hydraulic infrastructure. For instance, the present case shows that solution number #6 may be considered as one that satisfies almost simultaneously the three objective functions with quasi-equal weights of w 1 ¼ 0:32, w 2 ¼ 0:37 and w 3 ¼ 0:31. However, for examples, solution number #1 is the one that gives more importance to chute cost minimization with w 2 ¼ 0:59, but no importance to the stilling basin cost. Also, solution number #17 is the one where the energy dissipation (DE=E m ¼ 1 À f 1 ) within the chute is maximal, but with the largest spillway cost.
Comparatively to finding a unique optimal solution in single-objective optimization, multi-objective problems present a possibly uncountable set of solutions (Türkyılmaz et al. 2020). Providing a small number of representative optimal solutions may be of great help for project managers to make their decision. Retaining a solution to be realized depends on several parameters. For example, solution number #10 in Table 2 is identical to that found by Haddad et al. (2010) using HBMO, but it doesn't lead to a satisfactory energy dissipation rate (w 1 ¼ 0:07). The optimal design dimensions given by the reference are: b ¼ 10 m, h s ¼ 0:917 m , and L B ¼ 12:62 m. According to Table 2, this solution doesn't lead to a good energy dissipation rate and not even to the minimal spillway cost as in solution #2.
A compromise solution between energy dissipation and total cost may be selected from the set of clustered solutions. Figure 6(a) shows variations of the total cost versus energy dissipation as calculated from optimal solutions (blank circles). Some 13 other geometries outside the optimal Pareto set are also added to the plot in order to show that they necessarily lead to higher construction costs. In addition, Figure 6(a) shows also that the total construction cost C T ¼ f 2 þ f 3 is directly proportional to the desired energy dissipation rate: C T ¼ aDE=E m . A very good linear adjustment is obtained with a line slope   Figure 6(b) for e ¼ 0:5 m, e ¼ 0:75 m, e ¼ 1:25 m, and e ¼ 1:5 m. It was found that the slopes of the different adjustment lines may be expressed by: a ¼ 5:5e þ 2:8; this relies on the total cost to the energy dissipation rate as follows: Equation (11), valid for the case study considered herein, may be derived for a real spillway and provided to the project managers as a very practical tool.
Finally, five selected optimal solutions are analyzed using refined numerical modeling based on the finite volume method. These solutions are listed in Table 3 with the corresponding residual energy rates reported from the optimization procedure and that evaluated from the numerical simulations. The relative differences between optimal and numerical results are also given. They show that the results are close to each other except for solution #17, where the difference reaches 18%. This relatively great difference is due to the cluster center position of the solution. This center is situated near the transition flow zone just outside the border of the skimming flow zone. This does happen if no correction technique is used in the clustering phase to replace the centers at good positions. Solution #17 is thus outside the assumptions made at the beginning for application of Equations (2) and (4), so the right residual energy rate calculation is that given by the numerical simulation. Other numerical results related to field velocity, the shear stresses and flow turbulence may also be given if needed to enhance the decision making for the final design parameters. Therefore, the optimization procedure here is very useful; it leads to selection of very few optimal solutions to be analyzed with heavy numerical simulation or with expansive and time-consuming experimental modeling.

CONCLUSION
The most important points in hydraulic engineering dealing with spillways are reduction of both residual energy and total construction cost. In the present work, an optimization procedure is proposed for stepped spillway design that gives maximum energy dissipation rate and minimum construction cost considering independently the chute and stilling basin costs.

Corrected Proof
The multi-objective particle swarm optimization (MOPSO) is first employed to determine the best spillway dimensions (spillway width, step height and stilling basin length) that satisfy the above three objectives. A large number of Pareto solutions are obtained and then reduced to a small exploitable number by using K-means clustering algorithm. Then the clustered solutions are weighted following the pseudo-weight vector approach to facilitate the decision making. The obtained results are also used to derive an approximate expression for stepped spillways total cost versus flow energy dissipation rate. The optimization procedure may be considered as a useful tool to select a few optimal solutions to be analyzed with heavy numerical simulations or with experimental modeling for more desired results. It was found that, due to the implemented clustering technique, one among the solutions provided by the present procedure leads to an energy dissipation rate that differs from the one computed by numerical simulations. The clustering algorithm used herein needs enhancement to take into account flow regime limits in computing clusters means' positions. Moreover, including more flow characteristics in the optimization procedure, such as pressure distributions and shear stresses, will bring more information for the decision making about the optimal solution to be retained.
Finally, the present procedure may be applied to any multi-objective optimization problem, providing equations governing the objective functions to minimize. The developed code can be easily extended or modified to solve other optimization problems.

DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.