Application of genetic programming in shape optimization of concrete gravity dams by metaheuristics

A gravity dam maintains its stability against the external loads by its massive size. Hence, minimization of the weight of the dam can remarkably reduce the construction costs. In this paper, a procedure for finding optimal shape of concrete gravity dams with a computationally efficient approach is introduced. Genetic programming (GP) in conjunction with metaheuristics is used for this purpose. As a case study, shape optimization of the Bluestone dam is presented. Pseudo-dynamic analysis is carried out on a total number of 322 models in order to establish a database of the results. This database is then used to find appropriate relations based on GP for design criteria of the dam. This procedure eliminates the necessity of the time-consuming process of structural analyses in evolutionary optimization methods. The method is hybridized with three different metaheuristics, including particle swarm optimization, firefly algorithm (FA), and teaching–learning-based optimization, and a comparison is made. The results show that although all algorithms are very suitable, FA is slightly superior to other two algorithms in finding a lighter structure in less number of iterations. The proposed method reduces the weight of dam up to 14.6% with very low computational effort. Subjects: Computer Aided Design (CAD); Structural Engineering; Water Engineering


Introduction
Dams are structures built across a stream to form a reservoir. Dams are among the most important hydraulic structures used for various purposes such as water storage for domestic, industrial, and agricultural usage, flood controlling, generation of electricity, and so on. Gravity dams are solid concrete structures that maintain their stability against design loads from the geometric shape and the mass and strength of the concrete. Generally, they are constructed on a straight axis, but may be slightly curved or angled to accommodate the specific site conditions (USBR, 1976). Generally, the overall cost of dam construction is very high, and optimization methods are suitable tools for economically designing concrete gravity dams, owing to the fact that the shape of a gravity dam is so influential in its stability and that the overall cost of a gravity dam is proportional to the volume of concrete employed in dam construction. Many population-based optimization methods are recently used in engineering practice. Genetic algorithm (GA), particle swarm optimization (PSO), ant colony optimization and recently developed firefly algorithm (FA), and teaching-learning-based optimization (TLBO) are some examples. The main problem in utilizing such population-based approaches for solving real-life engineering problems is their high computational effort. First, a population of solution candidates (say 100) should be initialized and the solution is then sought iteratively (say for 100 iterations). Hence, the objective function should be evaluated several times (say 10,000). In engineering problems, the objective function is evaluated through analysis, e.g. by static or dynamic analysis of trusses, bridges, or dams which is itself time-consuming. This makes the overall procedure to be computationally expensive. A brief literature review on the application of three metaheuristics used in this study, i.e. PSO, FA, and TLBO, is presented in the next paragraph.
Optimal design of gravity dams is a fairly complicated procedure since many factors such as fluid-structure interaction, and seismic loads should be taken into account. Furthermore, some design constraints should be satisfied to ensure safety of dams against sliding and overturning as well as fracture. Several researches can be mentioned regarding optimization of concrete gravity and arch dams. Most of these researches are focused on design of optimal shape of arch dams rather than gravity dams. Simoes and Lapa (1994) and Simoes (1995) presented optimal design of dams using quadratic programming algorithm. They investigated the optimum shape of dams under static and dynamic loadings. The objective function was the weight of dam. They reported 13% reduction in the actual weight of dam after the optimization. Ghazanfari Hashemi, Bahraninejad, and Ahmadi (2009) employed simulated annealing (SA) approach for the shape optimization of gravity dams and compared their results with the results obtained by sequential quadratic programming (SQP) approach. They reported that the weight of dam can be reduced up to 30.18% using SA and up to 30.11% using SQP. Li, Jing, and Zhou (2010) used ANSYS software for the analysis of dam models and optimized a concrete gravity dam using genetic algorithm. They could reduce the weight of the dam up to 11.9% by this approach. Salmasi (2011) employed genetic algorithm in conjunction with the spreadsheet EXCEL to find the optimal shape of gravity dams. The static analysis was performed in this study. Qi (2012) used constrained nonlinear complex optimization algorithm to optimize the shape of gravity dams. He reported 16% reduction in cross-sectional area of the dam after optimization is accomplished.
In this study, shape optimization of concrete gravity dams is carried out using a rather different approach. The proposed method is a hybridization of metaheuristics and genetic programming (GP) to reduce the computational effort in the procedure of optimization. GP is a powerful tool which has been already employed to solve various problems in science and engineering (Danandeh Mehr, Kahya, & Yerdelen, 2014;Roushangar, Mouaze, & Shiri, 2014;Silva, Santos, Matos, & Costa, 2014;Valencia, Haak, Cotillon, & Jurdak, 2014). In order to introduce the approach used in this research, shape optimization of Bluestone dam is presented as a case study. The objective function is the weight of the dam, and design constraints are necessary safety factors against sliding and overturning and satisfying behavior constraint as well. The design variables consist of four geometric variables. The open access software CADAM is used for dam models. Despite being a very effective tool for the analysis of gravity dams, the software cannot be directly employed for optimization purposes since in its present form, it cannot be linked with programming languages such as MATLAB or FORTRAN. In order to overcome this problem, GP is employed. The proposed method has also reduced the computational effort required in metaheuristics. Several models of the dam including 322 models are analyzed using CADAM in which pseudo-dynamic analysis is performed for each model. The results of the analyses are used as a database to develop appropriate formulas for design constraints by means of GP. This procedure eliminates the requirement of structural analysis in every iteration in evolutionary optimization approaches. The method is combined with three metaheuristics including PSO, FA, and TLBO, and the results are compared.

Problem formulation
For stability requirements, the dam must be safe against overturning and sliding. Moreover, the safe stresses in the concrete of the dam or in the foundation material shall not be exceeded. The safety factor against overturning (OSF) is defined as the ratio between the resisting moments (M R ) and overturning moments (M O ): Many of the loads on the dam are horizontal or have horizontal components which are resisted by frictional or shearing forces along horizontal or nearly horizontal planes in the body of the dam, on the foundation or on horizontal or nearly horizontal seams in the foundation. A dam will fail in sliding at its base, or at any other level, if the horizontal forces causing sliding are more than the resistance available to it at that level. The resistance against sliding may be due to friction alone, or due to friction and shear strength of the joint. The stability of a dam against sliding is evaluated by comparing the minimum total available resistance to the total magnitude of the forces tending to induce sliding. Sliding resistance is a function of the cohesion inherent in the materials and at their contact and the angle of internal friction of the material at the surface of sliding. The factor of safety against sliding (SSF) can be computed from the following equation: ∑ F H is the sum of all the horizontal forces acting above the foundation level or at the level of the assumed sliding surface, is the internal friction angle, c is the cohesion, and A b is the area of the dam-foundation contact surface, or of the assumed sliding surface.
In shape optimization of gravity dams, the optimal shape of the dam with minimum weight should be determined. Hence, the objective function is the weight of the dam. Moreover, some design criteria as design constraints should be satisfied. In this study, optimization of dam is performed based on USBR standard. Based on USBR, the minimum required safety factor against sliding is 1.3 and the minimum required OSF is 1.1 for seismic loads. Furthermore, for satisfying behavior constraints, the principal stresses within dam body should not exceed the allowable tensile and compressive stresses of concrete (Beser, 2005). According to these design criteria, the optimization problem can be defined as follows: in which W is the total weight of the dam, c is the unit weight of concrete, V is the total volume of the dam, 1p and 2p are principal stresses, + and − are allowable tensile and compressive stresses, respectively, X is the vector of design variables, X L is the lower bound of the design variables, and X U is the upper bound. Since the value of c is constant, for unit width of the dam, the cross-sectional area of the dam can be minimized instead of W.
In this study, the allowable tensile and compressive stresses of concrete are considered to be less than the tensile and compressive strength of concrete as follows:

Structural analysis procedure for the case study
For analyzing gravity dams under hydrostatic and hydrodynamic pressures, gravity load, and seismic loads, various methods are available in the range of static to fully dynamic analysis. On the other hand, the structure should be analyzed many times either in the evolutionary optimization procedure or in artificial intelligence methods in order to find the optimal structure. Therefore, an effective and sufficiently accurate analysis approach should be chosen to reduce the computational effort. In this study, pseudo-dynamic analysis is selected for this purpose rather than complicated dynamic analysis. Pseudo-dynamic analysis has been found to be reliable in analyzing concrete gravity dams (Siamardi, 2007). This type of analysis is not only computationally more effective than dynamic analysis but also produces acceptable results close to those of dynamic analysis.
The pseudo-dynamic analysis is based on the simplified response spectra method. A pseudodynamic analysis is conceptually similar to a pseudo-static analysis except that it recognizes the dynamic amplification of the inertia forces along the height of the dam (Leclerc, Leger, & Tinawi, 2002). In this method, the dynamic magnification factor, inertial forces at the height of the dam, is considered. However, the oscillatory nature of the amplified inertia forces is not considered. That is the stress and stability analyses are performed with the inertia forces being continuously applied in the same direction. Since the pseudo-dynamic method does not recognize the oscillatory nature of earthquake loads, it is also appropriate to perform the safety evaluation in two phases: (1) the stress analysis using peak spectral acceleration values and (2) the stability analysis using sustained spectral acceleration values. It is assumed in these analyses that the dynamic amplification applies only to the horizontal rock acceleration. The period of vibration of the dam in the vertical direction is considered sufficiently small to neglect the amplification of vertical ground motions along the height of the dam (Leclerc et al., 2002).
The basic input data required to perform a pseudo-dynamic analysis, using the simplified response spectrum method proposed by Chopra (1998), are: (1) Peak ground and spectral accelerations, (2) Dam and foundation stiffness and damping properties, (3) Reservoir bottom damping properties and velocity of an impulsive pressure wave in water, and (4) Modal summation rule.
In this study, a free available code called CADAM is employed to perform pseudo-dynamic analysis of models. CADAM has been specifically developed for the analysis of gravity dams. This code was first presented by Leclerc et al. at Montreal University, Canada, in 2000(Leclerc et al., 2002Leclerc, Léger, & Tinawi, 2003). The main objective of the code is to analyze stability of concrete gravity dams against overturning, sliding, and breakage. It has been found that the code has great capability in analyzing gravity dams in spite of its simplicity. Because of its ability in analyzing gravity dams considering the most effective and influential forces, CADAM is one of the best codes in this area. CADAM analyzes dam models and directly reports the safety factors against sliding and overturning. It also reports the principal stresses developed within the dam body. Therefore, the implementation of the code is straightforward and no post-processing of the results is required after analysis.
The case study in the current article is devoted to optimizing the Bluestone Dam. This dam is located on the New River which is a tributary of Kanawha River 3.  Figure 1 along with other specifications reported in tables. Table 1 shows material properties considered in this study for dam and foundation; in Table 2, the properties of drainage gallery are given, and Table 3 reports the parameters used in seismic analysis of the dam. Figure 2 shows one of the models developed in CADAM, schematically.
In their study on the Bluestone Dam, Tekie and Ellingwood carried out a seismic hazard study on the dam, and the result is shown in Figure 3 (Takie & Ellingwood, 2003). Based on their study, the spectral acceleration of 0.36 g is considered related to a return period of 2,475 years in which g is the gravity acceleration. Moreover, based on a hydrologic analysis, if imminent failure flood occurs, the water depth in the reservoir will be 52 m (Takie & Ellingwood, 2002). In current study, it is assumed that the dam is full (water depth equals 48.15 m in the reservoir) and all analyses are performed based on this assumption.   ***Square-root-of-the-sum-of-squares of the first mode and static correction for higher modes.

The proposed method
The current study aims to hybridize the results found by GP with an evolutionary optimization method. For this purpose, dam models should be developed and the GP should be applied on the results in order to find appropriate relations for design constraints. Then, a suitable evolutionary algorithm can be used to determine the optimal values of design variables. The following sections illustrate this procedure in more detail.

Development of models and GP
A concrete gravity dam resists against the external loads by its own weight. For specified density of concrete, the weight of the dam is a function of its size. Hence, for other fixed material properties given in Table 1, the design criteria of dam including required safety factors against overturning and sliding and maximum principal stresses within the dam's body are related to the geometry of the dam. Figure 4 shows a schematic view of Bluestone Dam with geometric design variables considered in this study. As it is clear from the figure, maximum number of geometry variables allowed in CADAM (five geometric variables) is considered and optimized. The ranges of variables are considered as follows: In order to obtain a fairly complete database, 322 models of the dams with geometric values in the range of Equation 5 were developed and a complete analysis was carried out on each model by CADAM to obtain safety factors against overturning and sliding. Moreover, maximum principal stresses within dam body were obtained for all models. The results are then used to obtain appropriate formulas for these design parameters using GP. These formulas will be employed in optimization procedure in order to expedite the process of finding optimal values of design variables and to reduce computational cost of the approach. A glance on concepts of GP is presented in following paragraphs.
A branch of artificial intelligence is GP, which is an evolutionary algorithm-based methodology inspired by biological evolution to find computer programs that perform a user-defined task. In the 1960s and early 1970s, evolutionary algorithms became widely recognized as optimization methods. The first presentation of modern "tree-based" GP was introduced by Cramer (Takie & Ellingwood, 2003). This study was later greatly developed for application of GP in various complicated optimization and search problems. Recently GP has produced many novel and outstanding results in areas such as engineering, sorting, and searching, due to improvements in GP technology and the exponential growth in CPU power.
GP evolves computer programs, represented in memory as tree structures (Cramer, 1985). Trees can be easily evaluated in a recursive manner. Every tree node has an operator function and every terminal node has an operand, making mathematical expressions easy to evolve and evaluate  (see Figure 5). Therefore, GP favors the use of programming languages such as Lisp that naturally embodies tree structures. Table 4 gives the parameter setting for GP used in this study. As it is customary, 70% of data are randomly considered for training and the rest 30% for test.

Optimization methods
The proposed method can be effectively hybridized with any evolutionary algorithm. Three different metaheuristics including, PSO, FA, and TLBO are used to find the optimal shape of Bluestone Dam and the results are compared. In the following, brief reviews on algorithms are presented.

PSO algorithm
For the first time, Kennedy and Eberhart (1995) proposed the standard PSO based on the swarm behavior such as fish and bird schooling in nature. This approach searches a space of an objective function by individual agents called particles. A swarm is formed by the collection of particles. The particle movement has stochastic and deterministic components. The particle is attracted toward the position of the current global best while at the same time it has a tendency to move randomly. When a particle finds a location that is better than any previously found locations, it updates it as the new current best for particle i. There is a current best for all n particles. The aim is to find the global best among all the current bests until the objective no longer improves or after a certain number of iterations are carried out. In standard PSO algorithm, the swarm is updated by the following equations: Figure 5. Representation of a function as a tree structure in GP. where X k i and V k i represent the current position and the velocity of the ith particle at time k, respectively; P k i is the best previous position of the ith particle (called Pbest i ) at time k and P k g is the best global position among all the particles in the swarm (called gbest) at time k; r 1 and r 2 are two uniform random numbers between zero and 1.0; and c 1 and c 2 are two cognitive and social accelerating constants. Shi and Eberhart (1998) introduced an inertia term ω in Equation 6 and rewrote it to the form: They proposed that ω be selected such that 0.4 ≤ ≤ 0.9. Consequently, they report improved convergence rates when ω is decreased linearly during the optimization.
In this study, c 1 and c 2 were taken to be 1.0 and the inertia term begins with a value equal to 0.9 and it is decreased linearly to a value equal to 0.4 during the iterations.

Firefly algorithm
One of the latest metaheuristics algorithms is FA which is a nature-inspired algorithm developed by Yang (2009), inspired by the light attenuation over the distance and fireflies' mutual attraction. FA idealizes some of the characteristics of the firefly behavior in nature. They follow three rules: (1) all the fireflies are unisex, (2) attractiveness is proportional to their flashing brightness which decreases as the distance from the other firefly increases due to the fact that the air absorbs light.
Since the most attractive firefly is the brightest one, to which it convinces neighbors moving toward. In case of no brighter one, it freely moves any direction, and (3) brightness of every firefly determines its quality of solution; in most of the cases, it is proportional to the objective function.
During the loop of pair-wise comparison of light intensity, the firefly with lower light intensity will move toward the higher one. The moving distance depends on the attractiveness. The light intensity of each firefly is proportional to the quality of the solution, it is currently located at. In order to improve own solution, the firefly needs to advance towards the fireflies that have brighter light emission than is his own. Each firefly observes decreased light intensity, than the one the firefly actually emits, due to air absorption over the distance.
Attractiveness of a firefly is evaluated by the following equation (Yang, 2009): in which 0 is the attractiveness in distance r = 0 and is light absorption coefficient in the range [0, ∞). The distance r between firefly i and j at x i and x j is defined as the Cartesian distance: where x i,k is the kth component of the spatial coordinate x i of the ith firefly and d is the number of dimensions (Yang, 2009). Moreover, the movement of firefly i which is attracted by a more attractive or brighter firefly j is given by the following equation: in which α being the randomization parameter such that ∈ 0, 1 , and is a vector of random numbers drawn from a Gaussian distribution or uniform distribution in the range [0, 1]. For most problems 0 = 1 is considered. In the current study, the other parameters were set to = 0.1 and = 0.1.

TLBO algorithm
One of the latest metaheuristics is TLBO algorithm (Rao et al., 2011). Similar to most other evolutionary optimization methods, TLBO is a population-based algorithm inspired by the learning process in a classroom. The searching process consists of two phases, i.e. teacher phase and learner phase. In the teacher phase, learners first get knowledge from a teacher and then from classmates in the learner phase. In the entire population, the best solution is considered as the teacher (X teacher ).
On the other hand, learners learn from the teacher in the teacher phase. In this phase, the teacher tries to enhance the results of other individuals (X i ) by increasing the mean result of the classroom (X mean ) towards his/her position X teacher . In order to maintain stochastic features of the search, two randomly generated parameters r and T F are applied in update formula for the solution X i as: where r is a randomly selected number in the range of 0 and 1 and T F is a teaching factor which can be either 1 or 2: Moreover, X new and X i are the new and existing solutions of i, (Rao et al., 2011(Rao et al., , 2012. In the second phase, i.e. the learner phase, the learners attempt to increase their information by interacting with others. Therefore, an individual learns new knowledge if the other individuals have more knowledge than him/her. Throughout this phase, student X i interacts randomly with another student X j (i ≠ j) in order to improve his/her knowledge. In the case that X j is better than X i (i.e. f (X j ) < f (X i ) for minimization problems), X i is moved toward X j . Otherwise it is moved away from X j : If the new solution X new is better, it is accepted in the population. The algorithm will continue until the termination condition is met.

Constraint handling
The most common method for constraint handling in optimization algorithms is penalty function approach. The method has been already employed successfully to deal with constraints in other optimization problems (Ben Hadj-Alouane & Bean, 1997;Deb, 2000;Nanakorn & Meesomklin, 2001;Perez & Behdinan, 2007;Salajegheh & Khosravi, 2011). The main reason for the popularity of the method is its simplicity and its direct applicability regardless of the optimization method being used. This formulation utilizes general information of the collection of particles, such as the average of the objective function and the level of violation of each constraint in each iteration in order to define different penalties for different constraints as (Ben Hadj-Alouane & Bean, 1997): in which k i is the penalty parameter and it is calculated in each iteration by: with f (x) being the objective function and m being the number of constraints. Moreover, g i (x) is a specific constraint value so that violated constraints have values greater than zero, f (x) is the average of objective function in current particles and ḡ j (x) is the violation of the jth constraint averaged over the current collection of particles.
The illustration of Equation 16 is that the problem is actually solved as an unconstrained one, where in the minimization case, the objective function is designed such that non-feasible solutions are characterized by high-function values.

Results
The results of the implementation of aforementioned methodology including the results of GP implementation and the results of optimization are presented in following subsections.

Genetic programming
In the current study, the following operator functions and operands were used to find suitable expressions for safety factors against sliding and overturning according to the database developed by analyzing the models: After evolutions in GP, the following expressions for safety factor against sliding (SSF) and OSF, and principal tensile stress were found: Figure 6(a-c) shows the convergence history in GP for finding Equations 18,19,and 20, shows the effectiveness of Equations 18, 19, and 20 in predicting SSF, OSF and tensile stress respectively. Moreover, Table 5 reports some statistical criteria for testing goodness of the aforementioned equations in the prediction of design constraints. As Figure 7(a-c) as well as Table 5 show, the performances of Equations 18 and 19 in predicting SSF and OSF of the dam are excellent, whereas the performance of Equation 20 in predicting principal tensile stress is good. Therefore, these expressions are reliable to be used in the optimization procedure in the next step. It is worth pointing out that the analyses showed that the compressive principal stresses were found to be all well below the allowable stresses, and hence this criterion was not influential in the design of the current dam and no expressions for predicting this stress was required. (17)

Optimization
PSO, FA, and TLBO are used for finding the optimal values of design variables, whereas Equations 18 and 19 are used to evaluate safety factors against sliding and overturning, and Equation 20 for predicting principal tensile stress. Total numbers of 100 populations were used in each algorithm. Figure 8 compares the convergence history of minimizing cross-sectional area of dam using different metaheuristics, and Table 6 compares the results found by each algorithm. As Table 6 shows, all three methods lead to close results, indicating good performance of the algorithms. Optimal value of A = 892.8406 m 2 (W = 21,428.174 KN/m) was found by FA within 67 iterations, which shows slightly superior performance of FA compared to other two metaheuristics.
As Table 7 shows, for the design variables found by the algorithms, the values of SSF and OSF are the minimum values required for stability of the dam. After analyzing the optimal dam with CADAM, approximately the same values were found for SSF, OSF, and tensile stress, indicating accuracy of proposed formulas for these parameters. Moreover, the maximum tensile stress within the dam body was found to be 1p = 2,760 Kpa which is equal to the minimum allowable tensile stress; the maximum compressive stress is less than maximum allowable compressive stress within dam body, indicating stability of the dam against fracture as well. Compared to the actual weight of the dam   per unit length (W = 2,5101.44 KN/m), the optimal shape of dam shown in Figure 9 reduces the weight of dam up to a value equal to 14.6%. Figure 9 depicts the final shape of the cross section of the dam.

Conclusion
In this paper, the application of a powerful artificial intelligence approach, i.e. GP, in shape optimization of concrete gravity dam was introduced. A total number of 322 models of dams with regard to five design variables were developed and analyzed using pseudo-dynamic analysis. The results of the analyses were used to find appropriate formulas for design constraints. The accuracy of the formulas was verified, and they were used in an evolutionary optimization method in order to find optimal values of design variables. Three different metaheuristics, including PSO, FA, and TLBO were used for optimization and the results were compared. The procedure is computationally very effective since the necessity of time-consuming structural analysis in every evolution is eliminated. The results of applying the proposed method on a case study show that the procedure is capable of finding optimal values of design variables that satisfy all design constraints. One of the main aspects of this study is the applicability of artificial intelligence in solving real-life problems in engineering. In direct application of a structural analysis software to be used in a population-based optimization approach (say for 100 populations and 100 iterations), many analyses should be carried out to achieve the optimal design (say 10,000 analyses). Moreover, in some metaheuristics, the parameters of the algorithm should be tuned which means the overall procedure needs to be repeated several times. Hence, the populationbased algorithms are very cumbersome in dealing with huge structures or in the cases in which time-consuming dynamic analysis is necessary. In hybridization of artificial intelligence techniques with population-based optimization approaches, the number of such structural analyses is remarkably reduced.
In its current form, CADAM which was used in this study as the structural analysis platform does not permit more than two slopes to be defined in the downstream face of the gravity dam. In most cases, designers prefer to use three different slopes in the geometry of the dam at downstream face. This can be mentioned as the main shortcomings of this powerful software in the design of concrete gravity dams which may be resolved in the future updated versions.