Updated Kriging-Assisted Shape Optimization of a Gravity Dam

Abstract

The numerical simulation of the optimal design of gravity dams is computationally expensive. Therefore, a new optimization procedure is presented in this study to reduce the computational cost for determining the optimal shape of a gravity dam. Optimization was performed using a combination of the genetic algorithm (GA) and an updated Kriging surrogate model (UKSM). First, a Kriging surrogate model (KSM) was constructed with a small sample set. Second, the minimizing the predictor strategy was used to add samples in the region of interest to update the KSM in each updating cycle until the optimization process converged. Third, an existing gravity dam was used to demonstrate the effectiveness of the GA–UKSM. The solution obtained with the GA–UKSM was compared with that obtained using the GA–KSM. The results revealed that the GA–UKSM required only 7.53% of the total number of numerical simulations required by the GA–KSM to achieve similar optimization results. Thus, the GA–UKSM can significantly improve the computational efficiency. The method adopted in this study can be used as a reference for the optimization of the design of gravity dams.

1. Introduction

The shape of a dam is a major factor that influences its safety and economy. Therefore, designing a dam with a suitable shape is a significant problem in dam engineering [1]. Currently, the finite-element method is widely used in the optimization of gravity dam design [2,3,4]. However, a large number of numerical simulations are usually required when using optimization algorithms in finite-element software to optimize the shape of a gravity dam. To reduce the computational cost, some scholars have attempted to use surrogate models as alternatives to the finite-element analysis technique [5,6,7,8,9]

Among the numerous existing surrogate models, the Kriging surrogate model (KSM) [10] has been widely used because of its suitable performance. In the early stage of its development, the KSM served as an approximate alternative model to the time-consuming numerical analysis model for design optimization. Putra et al. [11] used the Kriging-based optimization method to optimize stent geometry design. Gaspar et al. [12] presented Kriging interpolation models as surrogate models for solving structural reliability problems involving nonlinear finite-element analysis. The suitability of the aforementioned optimization results depended largely on the accuracy of the KSM in the entire design space. In general, many training samples are required to increase the accuracy of a KSM.

A high sampling cost results in a limited number of samples being used for engineering optimization. A KSM constructed on the basis of a small number of samples may provide low-accuracy optimization results. In the updated KSM (UKSM), sampling is conducted only in the domain of interest. This strategy can considerably reduce the sampling cost while ensuring high prediction accuracy [13]. The UKSM has been successfully applied in many fields. Mohammadi-Amin et al. [14] used the KSM to estimate the aerodynamic coefficient trend for different spatial variables and effectively improved the accuracy of Kriging interpolation by adding accurate samples multiple times to the model. Augspurger et al. [15] constructed a UKSM response surface to develop practical thermal storage devices. Song et al. [16] presented an efficient optimization procedure for aerodynamic design. The selection of multiple samples according to infill criteria in each updating cycle resulted in the optimization procedure having suitable exploration ability. An efficient sequential surrogate-based optimization procedure was introduced in a previous study to optimize the design of Microridge Deep Drawing Punch [17]. However, the UKSM has been rarely used in the optimization of the shape of gravity dams.

In the present study, we adopted a structural optimization algorithm framework that couples GA and UKSM to optimize the shape of a gravity dam. Our main objective is to illustrate the efficient performance of GA-UKSM in the optimization of the gravity dam compared with GA-KSM. First, an optimization problem of a gravity dam was established. Second, the optimization problem was solved using a combination of GA and KSM (GA–KSM) and the GA–UKSM. The adopted structural optimization algorithm framework is expected to be a reference for the shape optimization of other dams.

2. Methodology

2.1. Formulaiton of the Optimization Problem

The cross-sectional geometry of the non-overflow section of a gravity dam is presented in Figure 1.

In Figure 1, A and H represent the dam crest width, dam bottom width, and dam height, respectively. Note that H₁ and H₂ are the heights of the upstream and downstream slope points of the dam, respectively. The parameters p (=H₁/H) and q (=H₂/H) represent the ratios of the heights of the upstream and downstream slope point vertices to the dam height, respectively. α₁ and α₂ are angels of upstream and downstream slope. The variables m₁ (=ctg(α₁)) and m₂ (=ctg(α₂)) represent the upstream and downstream slope coefficients, respectively. The objective of the gravity dam optimization problem is to minimize the cross-section area of the dam. The cross-section area can be express as follows:

$$S=\frac{{(pH)}^2{m}_{\mathit{1}}+{(qH)}^2{m}_{\mathit{2}}}{2}+AH,$$

(1)

where S is the cross-sectional area of the dam body; P, q, m₁, and m₂ are design variables that must be optimized; and the other symbols have the same meanings as previously mentioned. In this study, A and H are constants, so the objective function can be expressed as follows:

$$\min \frac{{(pH)}^2{m}_{\mathit{1}}+{(qH)}^2{m}_{\mathit{2}}}{2},$$

(2)

In the gravity dam optimization problem, behavior and stability constraints must be considered. These constraints are defined as follows:

$${\sigma }_{1,a}\le {\sigma }^{+},$$

(3)

$${\sigma }_{3,a}\le {\sigma }^{-} ,$$

(4)

$$k_a\ge K,$$

(5)

where σ_1,a and σ_3,a are the maximum and minimum principle stresses in the actual optimization process, respectively. For the entire earthquake process, the maximum value of σ_1,a and the minimum value of σ_3,a are considered. The parameters σ⁻ and σ⁺ represent the allowable compressive stress and allowable tensile stress, respectively; k_a is the actual safety factor against sliding; and K is the minimum required safety factor against sliding.

2.2. Solution of the Optimizaiton Problem

The GA–KSM and GA–UKSM were used to solve the optimization problem described in Section 2.1. The GA served as the optimization algorithm for solving the model. In the GA, individuals are evaluated by their fitness function values, and the fitness function values are obtained through numerical simulation. This process is computationally intensive; therefore, the KSM and UKSM were used to reduce the number of numerical simulations.

2.2.1. Genetic Algorithm (GA)

The GA is a heuristic algorithm that searches for optimal solutions by simulating natural evolutionary processes. It does not rely on gradients and has high robustness and global searching ability [18]. Therefore, the GA was employed to solve the constructed optimization problem. The design variables to be optimized were p, q, m₁, and m₂, which were encoded in binary in the GA. In this study, the penalty function [2] was used to satisfy the constraints involved in the optimization problem (Equations (3)–(5)). The fitness function containing the penalty term is expressed as follows:

$$F_{fit}=S+{10}^4\times [\left|\min \left\{k-K,0\right\}\right|+\max \left\{{\sigma }_1-{\sigma }^{+},0\right\}+\max \left\{{\sigma }_3-{\sigma }^{-},0\right\}],$$

(6)

where F_fit is the fitness function; k, σ₁, and σ₃ represent the safety factor against sliding, maximum principle stress, and minimum principle stress of each individual, respectively; and the other symbols have the same meanings as previously mentioned.

Multiple numerical simulations are usually required when the GA is executed using the ANSYS (Pittsburgh, Commonwealth of Pennsylvania, USA) software program for calculating the σ₁, σ₃, and k values of a gravity dam subjected to earthquake loading. To reduce the computational cost, the KSM was used in this study instead of numerical simulations.

2.2.2. Kriging Surrogate Model (KSM)

The KSM is a mathematical model for fitting discrete data. The essence of the KSM is to use the response information of known points to predict the response value of unknown points. The general steps for constructing the KSM are described in the following text. First, a sampling method (e.g., center composite design, orthogonal experimental design, uniform design, or Latin hypercube sampling (LHS)) is used to generate samples of design variables. Second, the samples are analyzed with a high-precision analysis model to obtain a set of input or output data. A fitting method is then used to fit the input or output relationship of these sample data to construct a surrogate model. Third, the fitting precision of the surrogate model is evaluated, and the output of the new design point is predicted.

• Sampling methods

Commonly used sampling methods include central composite design, orthogonal experimental design, uniform design, LHS, and inherited LHS (ILHS). LHS is the most widely used sampling method because this method is easy to implement and has satisfactory uniformity. The ILHS method [19] has the same uniformity as the LHS method but a higher sampling efficiency. Therefore, ILHS was used in the GA–KSM in this study. A total of 80% of the samples were used as training samples, and 20% of the samples were used as evaluation samples.

2.

Constructing the KSM

The general mathematical expression of the KSM is as follows:

$$\widehat{y}(X)=F(\beta ,X)+z(X),$$

(7)

where $\widehat{y}(X)$ expresses the deterministic response $y(X)\in I{R}^q$ for an n dimensional input $X\in D\subseteq I{R}^n$, F is the regression model which is a linear combination of p chosen functions f_j, and z(X) is a stochastic process.

$$\begin{array}{ll}F(\beta ,X)& ={\beta }_1{f}_1(X)+\dots {\beta }_p{f}_p(X)\\ & =[f_1(X)\dots f_p(X)]{\beta }_{:}\\ & =f{(X)}^T{\beta }_{:}\end{array}$$

(8)

The coefficients {β:} are regression parameters, and the random process z(X) is assumed to have mean zero, variance ${\sigma }^2$ and covariance between z(X_i) and z(X_j):

$$E[z(X)]=0,$$

(9)

$$Var\left[z(X)\right]={\sigma }^2,$$

(10)

$$Cov[z(X_i),z(X_j)]={\sigma }^{2}R(\theta ,X_i,X_j),$$

(11)

where R(θ, X_i, X_j) is a correlation function with parameters θ. The Gaussian function in Equation (12) is commonly used as a correlation function.

$$R(\theta ,X_i,X_j)={\displaystyle \prod _{k=1}^N{e}^{-{\theta }_k{\left|X_i^k-X_j^k\right|}^2}},$$

(12)

where X_i^k and X_j^k are the kth components of the training samples X_i and X_j, respectively, and θ_k is the constant of the correlation function in the kth direction of the sample. The linear weighted superposition interpolation of the response values of each sample is used to calculate the response value of the unknown sample as follows:

$$\widehat{y}(X)=w(X)Y,$$

(13)

where w(X) = (w₁, w₂, …, w_N)^T is the weighting factor to be determined and Y = (y₁, y₂, …, y_N)^T. According to the unbiased estimation condition, the following expression is obtained for the KSM through collation:

$$\widehat{y}(X)=f{(X)}^T{\beta }^{\ast }+r(X){\gamma }^{\ast },$$

(14)

where β* = (F^TR⁻¹F)⁻¹F^TR⁻¹Y and γ* = R⁻¹(Y − Fβ*). The parameters β* and γ* are fixed. The parameter r is expressed as follows: r = (R(θ, X, X₁), …, R(θ, X, X_N))^T. For more detailed derivation and application of the KSM, please refer to the reference [20].

3.

Accuracy evaluation

Accuracy evaluation is performed to judge whether a surrogate model is credible. If the accuracy meets the requirements, the surrogate model can replace the analysis model in the optimization process; otherwise, the surrogate model must be reconstructed. In this study, 20% of the samples were used to evaluate the KSM. The maximum absolute relative error (MARE) was used to assess the performance of the KSM. The MARE is defined as follows:

$$MARE=\max \left\{100\times \left|\frac{y_l-{\widehat{y}}_l}{y_l}\right|\right\}l=1, 2,\cdots N$$

(15)

where $y$ and $\widehat{y}$ are the actual value and predicted value, respectively, and N is the number of evaluation samples.

2.2.3. Updated Kriging Surrogate Model (UKSM)

The optimization results of the KSM largely depend on its approximation accuracy in the entire sampling space. To achieve the required approximation accuracy, a large number of samples must usually be collected; however, a large number of samples increases the calculation cost, which limits the application of the KSM in engineering optimization problems. By contrast, in the UKSM, in which a reasonable infill sampling criterion (e.g., minimizing the predictor (MP)) is constructed, a high approximation accuracy is only required near the optimal solution rather than in the entire design space. Thus, the UKSM can achieve a higher computational efficiency than the KSM.

In this study, MP [21] was used as the infill sampling criterion. The principle of the MP criterion is to find the minimum value of the objective function directly with a KSM and then take the point that minimizes the values of the objective as a new sample to update the KSM. Mathematically, the MP criterion can be expressed as follows:

$$\min \widehat{f}(X),$$

(16)

$$s.t.\{\begin{array}{l}{\widehat{h}}_n(X)\le 0 n=1,2,\dots N_c\hfill \\ X_L\le X\le X_U\hfill \end{array},$$

(17)

where X is a vector that contains the design variables p, q, m₁, and m₂; X_L and X_U are the lower and upper bounds of the design variables, respectively; $\widehat{f}\left(X\right)$ and $\widehat{h}\left(X\right)$ are a surrogate model of the objective function and a constraint function, respectively; and N_c is the number of constrained functions.

2.2.4. Optimization Procedure

• GA–KSM optimization procedure

To improve the sampling efficiency, samples were collected in the design space by using ILHS. Of the collected samples, 80% were used to obtain the ${\widehat{\sigma }}_1$, ${\widehat{\sigma }}_3$, and $\widehat{k}$ values of the KSM and 20% were used to evaluate the accuracy of the KSM. Subsequently, this study determined whether any of the MARE values between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, and k and $\widehat{k}$ were <1%. If none of the MARE values were <1%, 50 samples collected through ILHS were added to the sample set until the KSMs satisfied the accuracy requirements. The mathematical optimization model of the gravity dam was then solved using the GA. When convergence criterion (CC) ≤ 1‰ [Equation (18)], the optimization was terminated. A flowchart of the GA–KSM is illustrated in Figure 2a.

$$CC=\left|\frac{\frac{1}{5}{\displaystyle \sum _{gen=m}^{m+4}{S}_{gen}-\frac{1}{5}{\displaystyle \sum _{ge{n}^{\prime }=m+5}^{m+9}{S}_{ge{n}^{\prime }}}}}{\frac{1}{5}{\displaystyle \sum _{gen=m}^{m+4}{S}_{gen}}}\right|\times 100\% (m=1,2\dots ),$$

(18)

where gen and gen’ are the numbers of GA generations and S_gen and S_gen’ are the best objective values for gen and gen’, respectively.

2.

GA–UKSM optimization procedure

By using the MP criterion, new samples were added to the initial sample set to update the KSM [22,23,24]. First, the ${\widehat{\sigma }}_1$, ${\widehat{\sigma }}_3$, and $\widehat{k}$ values of the KSM were determined using 50 samples before optimization. Then, the GA was used to find the point minimizing the value of the objective function (when CC ≤ 1‰, the GA is considered to converge to the optimal solution), and the corresponding ${\widehat{\sigma }}_1$, ${\widehat{\sigma }}_3$, and $\widehat{k}$ values were obtained. Subsequently, the σ₁, σ₃, and k values corresponding to the optimal solution were calculated using ANSYS. If any relative error between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, or k and $\widehat{k}$ was >1%, the sample obtained through numerical simulation was added to the existing sample set to update the KSM until the entire optimization process converged. The schematic flowchart of the GA–UKSM is displayed in Figure 2b.

3. Case Study

3.1. Basic Information of the Gravity Dam

To investigate the computational efficiency of the GA–KSM and GA–UKSM for the shape optimization of a gravity dam, the Taolinkou Reservoir gravity dam was used as a test example. The basic section of the dam is displayed in Figure 1. The reservoir is located in the Qinglong River in the middle of the Yanshan settlement in Hebei, China. The maximum dam height is 89.24 m, the dam crest width is 8 m, the upstream water level is 84.74 m, and the downstream water level is 5 m.

In this study, four variables, namely p, q, m₁, and m₂, were defined for optimizing the gravity dam geometry.

Table 1. Lower and upper bounds of the design variables.

Design Variable	Lower Bound	Upper Bound
P	0.32	0.66
Q	0.86	0.93
m1	0.10	0.20
m2	0.60	0.80

presents the design variables of the optimization process that comply with China’s “Code for Design of Concrete Gravity Dams” and the engineering requirements.

In this study, we used the ground motion obtained from the record of an aftershock of magnitude 5.7 during the Tangshan earthquake that occurred in China on August 9, 1976. The seismic load included horizontal (Figure 3a and vertical Figure 3b) ground motions. The earthquake lasted for 5 s, and the seismic time step used in this study was 0.1 s. We ignored the quality of the foundation. To eliminate seismic wave propagation effects, we enter the seismic records on the surface of the dam foundation.

Under seismic loading, the allowable dynamic compressive strength of concrete was 30% higher than the allowable static compressive strength of concrete. The allowable dynamic tensile strength was 8% of the allowable dynamic compressive strength [25]. Moreover, the dynamic modulus of the dam body and dam foundation was 30% higher than the static elastic modulus. In line with China’s “Code for Seismic Design of Hydraulic Structures,” the minimum required safety factor against sliding (K) is 1.0. The mechanical properties of the dam body and foundation materials are presented in

Table 2. Material properties of the dam body and foundation.

Corresponding Area	Dam Body	Dam Foundation
Allowable value of static compressive strength (MPa)	11.17	14.91
Allowable value of dynamic compressive strength (σ−/MPa)	14.52	19.38
Allowable value of dynamic tensile strength (σ+/MPa)	1.16	1.55
Static modulus (GPa)	22.30	29.10
Dynamic elastic modulus (GPa)	29.00	37.83
Density (kg/m3)	2400	2500

.

3.2. Finite-Element Model of the Gravity Dam

In this study, finite elements were used to discretize the concrete and bedrock of the dam. Moreover, a two-dimensional finite-element model of the gravity dam (Figure 4) was established in ANSYS to calculate the maximum principle stress, the minimum principle stress, and the safety factor against sliding. The upstream and downstream of the dam foundation were expanded by 1.5 times the dam height, and the dam foundation was expanded downward by two times the dam height. In addition, normal constraints were applied upstream and downstream of the foundation, whereas fixed constraints were applied to the bottom of the foundation. The established model is based on the following assumptions: (1) the dam body and foundation are continuous, (2) the dam body and foundation materials are uniform, and (3) the concrete and bedrock of the dam body are isotropic elastic materials.

The loads acting on the dam are the gravity load, hydrostatic pressure, uplift pressure, and seismic load. The Westergaard additional mass method [26] was used to simulate the hydrodynamic pressure. The dynamic load effect was analyzed using the time-history analysis method [27,28,29]. Westergaard used the separation variable method and assumed that the dam surface facing the water is vertical, the dam body is rigid, and the reservoir water is a nonviscous liquid with small disturbances. He neglected the absorption of silt at the bottom of the reservoir, the surface microgravity wave of the reservoir water, and the compressibility of the reservoir water. Westergaard also assumed that the maximum hydrodynamic pressure on the dam surface has a parabolic relationship with the dam height. He obtained the hydrodynamic pressure of the calculation point using the following formula:

$$p=(7/8)a_x\rho \sqrt{hz},$$

(19)

where a_x is the horizontal acceleration, ρ is the water density, h is the designed water depth, and z is the distance from the calculation point to the water surface.

In the 1980s, Professor Clough [30] expanded the application scope of Westergaard’s additional mass method as follows:

$$M=(7/8)\rho A_i\sqrt{hz}{I}_i^T{I}_i,$$

(20)

where A_i is the subordinate area of the calculation point i on the dam surface, I_i is the normal vector of the calculation point, and the other symbols have the same meanings as previously mentioned.

3.3. Parameter Setting

In this study, the population size of the GA was set as 100 and the maximum generation as 500, and stochastic universal sampling was used as the selection operator. The single-point crossover probability was set as 0.7, and the mutation probability was set as 0.01. All the experiments were conducted on a computer equipped with an Intel Core i5-8250U 3.4 GHz CPU with 8 GB RAM (City of Santa Clara, Stage of California, USA).

3.4. Results and Discussion

When using the GA–KSM for optimizing the design of a gravity dam, the optimization results depend largely on the approximation accuracy of the KSM, which is directly related to the quantity of training samples. In this study, to meet the accuracy requirements of the KSM, the MARE threshold was set as 1%. As displayed in Figure 5, the MARE values between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, and k and $\widehat{k}$ decreased with an increase in the number of samples. When the MARE values were set to be ≤1%, the number of training samples was ≥1700. Such a high sampling cost is unacceptable in engineering optimization problems.

Therefore, we used the GA–UKSM to optimize the shape of the gravity dam. The relative errors between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, and k and $\widehat{k}$ varied with the number of iterations, as displayed in Figure 6.

As the number of iterations increased (i.e., as the number of new samples in the sample set increased), the relative errors between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, and k and $\widehat{k}$ decreased overall; however, volatility occurred. This volatility may have occurred because new optimal solutions were constantly being discovered during the process of adding samples to the sample set. When the optimal solution tended to be stable, the relative errors gradually decreased. This result indicates that the application of the MP criterion causes the UKSM to not have a high approximation accuracy in the entire design space. The model has a high approximation accuracy only near the optimal solution. Finally, only 78 sample points must be added to the KSM to meet the convergence condition.

Table 3. Comparison of the optimization results obtained with the GA–KSM and GA–UKSM.

Optimization Methods	Design Variable				H1(m)	H2(m)	S(m2)	Number of Simulations
	p	q	m1	m2
GA–KSM	0.395	0.872	0.195	0.682	35.25	77.82	2899.4	1700
GA–UKSM	0.396	0.872	0.186	0.680	35.34	77.82	2888.9	128

presents a summary of the optimization results obtained for gravity dam design with the GA–KSM and GA–UKSM.

As presented in Table 3, the difference in the cross-sectional area of the dam body obtained with the two algorithms was approximately 0.36%. In theory, both the compared algorithms can be used for the optimization of gravity dam design. However, when using the GA–KSM to optimize the gravity dam shape, at least 1700 samples are required in the KSM sample set to make the MARE values between σ₁ and ${\widehat{\sigma }}_1$, σ₃ and ${\widehat{\sigma }}_3$, and k and $\widehat{k}$ of <1%. Thus, at least 1700 numerical simulations are required in the GA–KSM (Figure 5). Such a high sampling cost is unfavorable in gravity dam optimization design. By contrast, the number of numerical simulations required for the GA–UKSM (128 numerical simulations, i.e., 128 samples, including 50 samples required to construct the initial KSM and 78 samples obtained using the GA to improve the accuracy of the KSM according to the MP criterion) was only 7.53% that required for the GA–KSM. Thus, the GA–UKSM considerably outperformed the GA–KSM in terms of the optimization efficiency.

Moreover, the differences in the design variables corresponding to the optimal solutions of the two algorithms were not obvious, which indicates that the UKSM can guarantee higher precision with fewer samples than the KSM can. The aforementioned finding also indicated that the optimization process of the GA–UKSM converged to a stable solution (Table 3).

In this study, only one sample that did not satisfy the convergence condition was added to the sample set to update the KSM in an updating cycle. In future research, an algorithm that simultaneously adds N (N > 1) samples to the sample set should be considered to improve the probability of obtaining an optimal solution. In addition, three-dimensional problems where spanwise dam shape and curvature and changes of dam foundation elevations in both axial and spanwise directions can be taken into account.

4. Conclusions

A combination of the GA–UKSM was adopted to optimize the shape of a gravity dam subjected to earthquake loading. To improve the optimization efficiency and prediction accuracy of the KSM, the MP criterion was used to update the KSM in each updating cycle. The optimization efficiency of the GA–UKSM was approximately one order of magnitude higher than that of the GA–KSM. As the samples in the sample set increased, the prediction accuracy of the UKSM increased. Consequently, the GA–UKSM outperformed the GA–KSM in terms of the computational efficiency and number of numerical simulations. The optimization design method for gravity dam adopted in this research can provide indirect guidance for engineering designers in hydraulic structure design, improving their optimization design efficiency.