Multi-parameter inversion of concrete face rockfill dam using wild horse optimizer and optimal polynomial chaos kriging

2.1. Optimal polynomial chaos kriging

The basic idea of PCE is to represent the output response as a linear combination of polynomial basis functions. The input variable $X$ is a random vector with independent components, and its distribution is governed by the joint Probability Density Function (PDF) $f_X$. The PCE is defined as:

where $y_{\alpha }$ represents the expansion coefficients, ${\mathrm{\Psi }}_{\alpha }\left(X\right)$ are multivariate orthogonal polynomials.

To simplify the computation, a truncated PCE is used [36, 37]:

where $A\subset N^M$ is the set of truncated polynomials.

The error ${\epsilon }_p$ is estimated using the leave-one-out (LOO) cross-validation [38, 39]:

where ${\mu }_{{\widehat{Y}}^{\left(PCE\right)},(-i)}$ represents the mean of the PCE model values obtained using all experimental design points $X$, excluding $X_i$.

Kriging is an interpolation method based on Gaussian processes, where the output is a sum of a trend term and a variance term [40]:

where the ${\beta }^{T}f\left(x\right)$ represents the mean value (trend) of the Gaussian process, which is a linear combination of polynomial basis functions $f\left(x\right)$ and regression coefficients $\beta$. ${\sigma }^{2}Z(x,\omega )$ is a Gaussian process that models the correlation between input samples, it includes the variance ${\sigma }^2$ of the Gaussian process and a zero-mean stationary Gaussian process $Z(x,\omega )$. The optimal hyperparameters are determined via Maximum Likelihood Estimation (MLE) [17] and used to predict responses for new sample points.

OPCK combines the global approximation of PCE and the local approximation of Kriging for a more accurate surrogate model. OPCK is defined as a universal Kriging model, with the trend consisting of a set of standard orthogonal polynomials [41-43]:

where ${\sum }_{\alpha \in A}{y}_{\alpha }{\mathrm{\Psi }}_{\alpha }\left(x\right)$ represents the weighted sum of standard orthogonal polynomials that describe the trend of the PCK model, while ${\sigma }^{2}Z(x,\omega )$ represent the variance of the Gaussian process as described in Section 2.2.

In OPCK, the optimal polynomial set is determined using the Least Angle Regression (LAR) process [44] within the PCE framework [45]. Trend and correlation parameters are computed using Kriging equations [46, 47]. The LAR algorithm ranks the polynomials based on their relevance to the residuals at each iteration. Finally, different PCK models are compared using the LOO error (Eq. (3)), and the model with the smallest LOO error is selected. A flowchart of the OPCK algorithm is provided by Schobi, Sudret, and Wiart [41].

2.2. Wild horse optimizer

This section describes the WHO model, inspired by wild horses’ social behaviors, as proposed by Naruei in 2021 [48]. WHO mimics grazing, chasing, dominance, leadership, and mating behaviors.

(1) Creating the Initial Population. The initial population of size N is divided into G groups, with PS representing the percentage of stallions (leaders). The remaining N−G members (mares and foals) are distributed among these groups.

(2) Grazing behavior. Mares and foals spend time grazing near the group leader. Their position is updated as follows:

where $X_{i,G}^j$ is the current position of the group members, $\mathrm{S}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{o}{\text{n}}^j$ is the position of the stallion (group leader), $Z$ is the adaptive mechanism. $R$ is a random number uniformly distributed within [-2, 2]. The grazing radius varies due to the cosine function.

(3) Mating behavior. Foals reach puberty and leave their groups to form or join new ones for mating. This behavior is modeled as:

the position of the new horse $X_{G,K}^p$ depends on the mating of group $i$ and group $j$ horses.

(4) Leadership behavior. The stallion leads the group to new areas, and the rest avoid it. This is modeled as:

where $\overline{\text{Stallio}{\text{n}}_{G_i}}$ is the new location of the group $i$ stallion, WH is the current location of the explored area, and $\mathrm{S}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{o}{\text{n}}_{G_i}$ is the current location of the group $i$ stallion.

(5) Leader Exchange and selection. The stallion is initially random, and as iterations proceed, the leader’s position is determined by fitness. If a member is better adapted, their positions are swapped:

The searching individual (the stallion) is positioned close to the optimal position as it is iteratively replaced.

3.1. Framework of method

The proposed WHO_OPCK method for parameter inversion is structured in four phases. Initially, a FE model is established using benchmark parameters to simulate the settlement behavior of a real structure, accounting for uncertainties through a probabilistic model. Following this, a sensitivity analysis is conducted using a surrogate model to identify critical parameters. The surrogate model is then combined with the WHO algorithm to minimize an objective function, ultimately refining the parameter inversion for accurate model updates. The framework of the WHO_OPCK method is illustrated in Fig. 1.

3.2. Procedure of WHO_OPCK

The procedure of implementing the WHO_OPCK method is detailed below.

Phase One – Define the problem (Step A): The first step is to establish a FE model using benchmark parameters to characterize the settlement behavior of the real structure. It is necessary to define a prior probabilistic model for the unknown parameters to be determined inversely, including distribution type, mean, and deviation. Since real structure parameters are often random, incorporating uncertainty in the input parameters is more realistic.

Phase Two – Parameter sensitivity analysis based on surrogate modeling (Step B~C): Create a few input parameter datasets randomly using Latin Hypercube Sampling (LHS) method. Perform Finite Element Analysis (FEA) and extract the corresponding output response dataset. Construct the OPCK model from the “input-output” data and assess its accuracy with Leave-One-Out (LOO) cross-validation error. If the accuracy is satisfied, perform a parameter sensitivity analysis based on the Sobol’ indices and the constructed OPCK model to filter out performance parameters to be inverted; If the accuracy is insufficient, generate additional sampling points.

Phase Three – Coupling the OPCK model with WHO and minimizing the objective function (Step D~F): According to the results of the sensitivity analysis, a certain number of inversion parameter datasets are randomly generated using LHS method. Subsequently, the FEA is executed, and the output response dataset is extracted. Then, the same method as in step C is used to construct the OPCK model and evaluate its accuracy. Once an accurate OPCK model is established, it can replace the FE model for calculating structural responses, bypassing the high computational cost of FE model. Finally, the OPCK model is coupled with the WHO algorithm by integrating the predicted and target (virtual measured) outputs into an objective function.

Phase Four – Parameter inversion post-validation (Step G): The optimal values from the parameter inversion are processed and used as actual input parameters to update the FE model.

The symmetric mean absolute percentage error (SMAPE) is a common metric employed in statistics for evaluating prediction accuracy. Unlike Root Mean Square Error (RMSE), it is expressed as absolute percentages, making it easier to interpret. In this paper, SMAPE is chosen as the objective function to access the accuracy of the WHO_OPCK method. The formula is as follows:

where $Y_{vm,i}$ is the virtual measured value, $Y_{OPCK,i}$ is the predicted value from the OPCK model, and $N$ is the number of fitted points.

4.1. Description of ZJTT CFRD

The ZJTT CFRD is the dam of the lower reservoir of a pumped storage power station. As illustrated in Fig. 2(a), it has a crest elevation of 209 m, a base elevation of 138 m, a maximum height of 71 m, and a total crest length of 404.25 m. The thickness of the concrete face slab is given by $T=$0.4 + 0.003$H$, where $H$ is the height from the calculated section to the top of the face slab, with a slope of 1:1.405. This study focuses on the period immediately following the completion of dam construction. During this period, the dam’s loading is predominantly influenced by its self-weight. As shown in Fig. 2(b-d), the displacement of the measurement point is used as the QoIs.

4.2. Simulation process

The layered stacking process of the rockfill dam was simulated using Abaqus, incorporating the Duncan E-B UMAT subroutine, developed in Fortran [21], for the rockfill material. Goodman contact elements with zero thickness were used to model the contact between the concrete face slabs and the bedding layer [28]. A three-dimensional numerical model was constructed to determine the displacement and stress distribution patterns of the dam body during the layered construction of the rockfill, as well as the stress-strain distribution patterns of the concrete face slab during the same period.

Fig. 1Framework to describe the WHO_OPCK method

The material zoning consists of five parts: main rockfill, secondary rockfill, transition, cushion, and concrete face slab, as shown in Fig. 3(a). The main rockfill, secondary rockfill, transition and cushion layers were modeled using Duncan E-B model with the material parameters listed in Table 1. The concrete face slab was modeled with linear elasticity and the material parameters are shown in Table 2. The Goodman contact element parameters in Table 3 reflect the friction and stiffness of the contact surface, ensuring reliable simulation results.

The stacking process was simulated using Abaqus life and death units, with the rock mass stacking divided into 15 steps, one step each for the cushion layer and the concrete face slab, as shown in Fig. 3(b).

Fig. 2Dam size and measurement profiles: a) longitudinal view of the dam, showing the profile locations where measurement points are distributed; b) cross-section at y= 361.99 m, showing the measurement point for Quantity of Interest QoI 1; c) cross-section at y= 229.99 m, showing the measurement points for QoIs 2, 3, 4, and 5; d) cross-section at y= 88.99 m, showing the measurement point for QoI 6

a) Longitudinal view of the dam

b) Cross-section at $y=$ 361.99 m

c) Cross-section at $y=$ 229.99 m

d) Cross-section at $y=$ 88.99 m

Fig. 3ZJTT CFRD: a) Material zoning layout of the dam, showing the distribution and classification of different materials used in the construction; b) Layered construction process of the dam, illustrating the sequence and method by which the dam is built layer by layer

a) Material zoning

b) Layered construction process

Given that the dam is still under construction, the measured data could not be obtained. To demonstrate the efficiency of the proposed WHO_OPCK method in engineering applications, a “hypothesis testing” strategy was employed to verify the method’s feasibility. First, the LHS method was used to randomly generate 20 sets of material parameters for the CFRD model. The corresponding displacements at the measurement points were then calculated and averaged based on FEA. To account for the discrepancies between numerical simulations and real-world measurements, 20 % Gaussian noise was introduced into the mean values. Finally, one set of noise-containing data was randomly selected as the virtual measurement data, as described in Step F in Fig. 1.

5.1. Parameter sensitivity analysis

5.1.2. Construction of the surrogate model

This section explores how the number of experimental designs ($N_{DOE}$) affects model accuracy and compares the OPCK surrogate model with classical PCE and Kriging. Surrogate models were built using varying $N_{DOE}$ to determine the optimal values for the OPCK model. A total of 140 initial samples were taken for the unknown parameters using LHS. Corresponding response output (QoI 1~6) was calculated. Then, surrogate models with different sizes of $N_{DOE}$ (i.e., 10, 20, 30, 40, 50, 60, 70, and 80) were constructed, this yielded 24 surrogate models. The $Er{r}_{LOO}$ for the six output responses in each case are shown in Fig. 5.

Fig. 4Seven-parameter random sampling results

From Fig. 5 the following results can be derived:

1) All three surrogate models show high accuracy across varying $N_{DOE}$. ut Kriging has lower accuracy and greater discrepancies compared to PCE and OPCK. The OPCK model outperforms PCE in accuracy.

2) For the PCE model, $Er{r}_{LOO}$ decreases significantly from $N_{DOE}=$10 to 40, then plateaus, showing no further improvement. In the OPCK model, the reduction is more pronounced between $N_{DOE}=$10 and 40, after which accuracy stabilizes for $N_{DOE}>$ 40.

3) While Kriging’s $Er{r}_{LOO}$ decreases with increasing $N_{DOE}$, it doesn’t reach the desired accuracy (below 5 %) until $N_{DOE}=$ 80.

Fig. 5The impact of the number of experimental designs (NDOE) on output accuracy and a comparison of three surrogate models

5.1.3. Sensitivity analysis

The Sobol’ index method, based on variance decomposition, determines the sensitivity of input variables by calculating their contribution to the total output variance [49]. In this study, the ErrLOO of the OPCK model is less than 10-2 when $N_{DOE}=$ 50, meeting the predictive accuracy requirements for large hydraulic structures. Hence, parameter sensitivity analyses were conducted using the OPCK model instead of the computationally expensive ZJTT CFRD model. The sensitivity analysis results, considering the impact of the seven unknown input parameters on the ZJTT CFRD model’s output, are shown in Fig. 6. Key observations include:

1) $\phi$, $K_b$, $K$, and $R_f$ exhibit high sensitivity in various QoI metrics, with $\phi$ and $K_b$ being particularly influential. These parameters are crucial in CFRD deformation and should be prioritized for calibration and optimization.

2) In contrast, $n$, $\mathrm{\Delta }\phi$, and $m$ show low sensitivity and weak effects on QoI metrics, contributing minimally to CFRD deformation and are of secondary importance.

3) Sensitivities of QoI indicators to the same parameter vary. For example, $\phi$ is significant in QoI 4 and 5, while $K_b$ is prominent in QoI 6, reflecting the system’s nonlinear response. Multidimensional indicators should be considered in parameter analysis.

Fig. 6Results of parameter sensitivity analysis based on the OPCK model

a) First-order Sobol’ indices

b) Total Sobol’ indices

The sensitivity analysis revealed that $\phi$, $K_b$, and $R_f$ are crucial for CFRD. The next step is multi-parameter inversion for these parameters to improve model accuracy and prediction.

5.2. Multi-parameter inversion

After the sensitivity analysis, $\phi$, $K_b$, $K$, and $R_f$ were identified as key parameters. We will now perform a multi-parameter inversion to calibrate the model and match the virtual measured data.

Section 5.1.2 discussed key factors influencing surrogate model accuracy. This section focuses on the impact of the number of search individuals ($N$) in the WHO_OPCK method, with Niter set to 500. Results are shown in Fig. 7, and observations are as follows:

1) Fig. 7(a) shows that as $N$ increases, convergence speed improves, particularly for $N\ge$50, where the error decreases rapidly. Smaller $N$ values (e.g., $N=$30 or $N=$35) lead to slower convergence, especially at the start, due to limited search space exploration.

2) The SMAPE value shows significant variation at lower $N$, dropping sharply at $N=$35 and $N=$40. For $N\ge$40, SMAPE stabilizes, indicating improved convergence accuracy as more individuals help find the global optimum.

3) Fig. 7(b) shows that computation time increases with $N$, as more individuals require processing more data. The increase stabilizes at $N=$60, where time overhead is moderate and accuracy remains high, offering an optimal balance.

Based on the results, surrogate models and inversion methods with $N_{DOE}=$ 50 and $N=$ 60 were chosen for comparison.

For the computationally expensive 3D CFRD, using $N_{DOE}=$ 50 for the surrogate model results in a total computation time of about 28 hours (0.56 hours per run) in ABAQUS. In contrast, the traditional WHO algorithm with direct FE model calls for iterative inversion would take approximately 16,800 hours (0.56×60×500) for $N=$60 and $N_{iter}=$500. This immense computational demand makes the traditional method impractical for 3D CFRD parameter inversion, limiting its real-world engineering application.

Fig. 7The effect of different N on the accuracy of WHO_OPCK method

a) Convergence properties

b) Computational time and convergence accuracy

The WHO optimization algorithm was combined with various surrogate models to compare the convergence of WHO_OPCK with WHO_SPCK, WHO_PCE, and WHO_Kriging. Fig. 8 shows the convergence speeds, accuracies, and computation times, while Table 5 presents the inversion results. The following conclusions can be drawn:

1) WHO_OPCK achieves faster convergence and greater robustness in parameter inversion. It outperforms other methods, especially WHO_SPCK and WHO_Kriging, by finding better solutions in fewer iterations while maintaining stability.

2) WHO_OPCK provides more accurate inversion results than WHO_SPCK, WHO_PCE, and WHO_Kriging, demonstrating superior convergence accuracy.

3) WHO_OPCK, though slightly slower than WHO_PCE, outperforms WHO_SPCK and WHO_Kriging in efficiency, requiring only 46 % and 16.8 % of their respective computation times, achieving an optimal balance between time and accuracy.

By combining the strengths of PCE and Kriging, WHO_OPCK excels in both accuracy and speed while keeping computation time reasonable, offering better cost-effectiveness in practical applications.

Meanwhile, we also coupled OPCK with classical optimization algorithms (PSO, SMA, SSA, WOA, GWO) [11, 12, 50-52] to compare convergence and computational efficiency. Fig. 9 shows the convergence speed, accuracy, and computation time, while Table 6 presents the inversion parameter results. The following conclusions can be drawn:

1) WHO_OPCK outperforms OPCK with classical algorithms in both convergence speed and accuracy, showing strong performance in complex CFRD optimization.

2) Fig. 9(b) shows that parameter inversion accuracy and convergence speed do not directly correlate with inversion time. Despite its high computation load, SSA does not achieve the highest accuracy, while WOA, with the second-highest time, has the lowest accuracy.

3) WHO_OPCK achieves the highest accuracy and shortest inversion time, making it the most efficient and effective method among the comparisons.

Fig. 8Comparison of WHO combined with different surrogate models

a) Convergence properties and accuracy

b) Computation time

Fig. 9Comparison of different optimization algorithms combined with OPCK

a) Convergence characteristics and accuracy

b) Computation time

The parameter inversion results are used to update the FE model, which is then simulated with the new parameters. The output is compared to the virtual measured values, with errors assessed using RMSE and SMAPE. The detailed output and error results are shown in Table 7 and 8, and Fig. 10 provides a visual comparison. The following conclusions can be drawn:

1) Fig. 10(a) shows that the FE model updated with WHO_OPCK achieves higher accuracy, with significant advantages in both RMSE and SMAPE. Among the surrogate models, WHO_Kriging has the greatest error and lowest accuracy.

2) Fig. 10(b) shows that the FE model updated with WHO_OPCK has better accuracy. SMA and WHO show similar errors, while WOA performs the worst, with the largest error and lowest accuracy.

Fig. 10Comparison of FEA errors based on inversion results and virtual measured data

a) WHO combined with different surrogate models

b) Different optimization algorithms combined with OPCK