An intelligent reliability-based design optimization of complex engineering structures using manifold learning-enhanced surrogate modeling

Abstract

With the increasing integration of computer science and engineering optimization, the capabilities of advanced machine learning in data processing and surrogate modeling for complex simulations have provided a transformative development direction for enhancing the computational performance of reliability-based design optimization (RBDO). However, when addressing highly nonlinear and complex problems, machine learning models often encounter challenges arising from the trade-off between modeling accuracy and computational efficiency during surrogate simulation analysis. Particularly in high-dimensional physical simulations, traditional machine learning approaches struggle to construct precise surrogate models, thereby limiting their effectiveness within the RBDO framework. Manifold learning represents an effective dimensionality reduction technique. To address these challenges, this study proposes a manifold learning-integrated surrogate modeling methodology along with a novel RBDO framework. The approach employs active dimensionality reduction to enhance the adaptability of advanced surrogate models across diverse complex engineering structures, complemented by sophisticated metaheuristic algorithms to improve optimization accuracy and efficiency. The proposed reliability-based design optimization method consists of three sequential stages: (1) developing high-fidelity surrogate models through manifold learning to replace complex engineering simulations; (2) conducting inverse reliability analysis incorporating the advanced surrogate model to transform uncertainty optimization problems into deterministic ones; and (3) solving for optimal design solutions using state-of-the-art metaheuristic algorithms. Experimental results demonstrate that the proposed RBDO framework reduces computational time by ∼68% compared to conventional design methods, validating the effectiveness, and superiority of this optimization approach for addressing complex engineering problems.

Introduction

The ever-increasing demands in contemporary society have led to products becoming increasingly complex, posing new challenges for machine learning-based reliability optimization design methods. In the engineering design of complex structures, various uncertain factors often come into play, such as material property variations, manufacturing process inconsistencies, and load fluctuations.^1,2 The large number of uncertainty factors may result in high costs for surrogate model construction and reduced accuracy of design outcomes. ³ Therefore, enhancing the precision and efficiency of surrogate model construction holds significant importance for uncertainty optimization design.

From the perspective of the theoretical framework for reliability design optimization, its methods can be categorized into two major types: time-dependent dynamic reliability optimization and time-independent static reliability optimization. ⁴ The former primarily focuses on the reliability characteristics of systems or products that evolve over time during their service life, while the latter concentrates on optimizing static variables such as structural parameters, geometric dimensions, and material properties during the design phase. The aim is to ensure, through deterministic or probabilistic methods, that the reliability metrics of products meet the requirements under given operating conditions. Currently, although dynamic reliability optimization holds significant value in long-term service scenarios, its computational cost and data requirements are substantially higher than those of static methods due to the need to process high-dimensional time-varying data and construct complex dynamic models. In contrast, time-independent reliability design optimization methods, being directly linked to the initial configuration of product design parameters, can more efficiently address the issues of uncertainty propagation and accumulation in complex structures. Moreover, these methods are more readily integrated with existing CAD/CAE tools, making them a focal point of research in both academia and industry. Influenced by uncertainties, structural reliability theory has played a pivotal role in advancing optimization design methodologies. Reliability-based design optimization (RBDO) has evolved from deterministic design methods by integrating mathematical programming techniques with reliability analysis theories.^5,6 RBDO is fundamentally a double-loop nested optimization issue. In this context, the outer loop is dedicated to the optimization of design variables, while the inner loop is responsible for conducting reliability analysis. Among the methods for tackling RBDO problems, the nested optimization approach stands out as the most straightforward. Within this approach, the reliability index approach (RIA)^7,8 and the performance measure approach (PMA)^9,10 are the two most frequently employed nested optimization techniques. To enhance computational efficiency, researchers have also proposed decoupling methods and single-loop methods for RBDO. The core concept of the decoupling method lies in approximating the failure probability function or the reliability index function by means of an inexpensive surrogate model prior to design optimization. Subsequently, the probabilistic constraints are substituted with this surrogate model, effectively converting the RBDO problem into an equivalent deterministic optimization problem. Utilizing the Karush–Kuhn–Tucker optimality conditions, the single-loop method^11,12 substitutes the probabilistic constraints, thereby converting the nested optimization problem into a single-loop optimization problem. This approach is particularly effective for problems with low nonlinearity in the limit-state functions. In recent years, more advanced optimization design methods have been developed based on traditional approaches. Chen et al. ¹³ proposed a decoupled optimization design framework incorporating a reliability optimization integration method, where the integral-based reliability analysis approach significantly improves the accuracy and efficiency of optimization design. Miska and Balzani ¹⁴ introduced a design optimization method under polymorphic uncertainties by integrating an extended optimal uncertainty quantification framework into the double-loop RBDO, thereby enhancing the optimization capabilities for complex engineering structures. Zhang et al. ¹⁵ proposed an RBDO method based on transfer learning, which transforms repetitive reliability analyses at different design points into a series of target-domain tasks, leading to substantial improvements in computational efficiency. Zhao et al. ¹⁶ presented a reliability optimization design method that combines probabilistic uncertainty and reliability analysis, successfully applying it to the optimization design of rooftop photovoltaic systems.

During the optimization process, as product complexity increases, evaluating limit states using numerical simulations or emulation methods often incurs significant time costs. To address this, the use of machine learning methods to construct surrogate models for auxiliary optimization design has emerged as a research focus.^17–19 By mapping input variables to corresponding output variables through mathematical expressions and establishing a mathematical relationship between them, optimization time and cost can be substantially reduced. Commonly employed surrogate models include the response surface method, ²⁰ polynomial chaos expansions (PCE), ²¹ support vector machine, ²² and Kriging,^23,24 among others. Pang et al. ²⁵ proposed a Kriging-based RBDO framework that balances exploration of the objective space by integrating global and local optimization modules, thereby enhancing optimization efficiency. Shi et al. ²⁶ introduced a scalable robust design optimization framework that constructs a Kriging model for system reliability indices through a two-stage active learning strategy, offering a novel approach for design optimization. Gargama et al. ²⁷ combined artificial neural networks with genetic optimization algorithms to optimize electromagnetic shielding structures, demonstrating that this method effectively improves structural reliability. Yu et al. ²⁸ proposed a structural optimization design method based on radial basis function neural networks and applied it to optimize offshore wind turbine support structures, showing significant improvements in optimization efficiency. Hamza et al. ²⁹ presented an RBDO method for mechanical engineering problems by integrating reliability design space techniques with an efficient hybrid algorithm combining adaptive differential evolution and Nelder-Mead local search, revealing that this approach effectively reduces the computational cost of RBDO. Yang et al. ³⁰ proposed a single-loop reliability-based design optimization (RBDO) framework based on a hybrid adaptive Kriging method. This approach effectively addresses the issues of computational inefficiency and convergence difficulties encountered in traditional single-loop methods. Additionally, they introduced ³¹ a single-loop RBDO framework that integrates an efficient single-loop local adaptive Kriging approximation, further enhancing (or more precisely, “reducing the computational cost for”—note the original Chinese has a potential logic issue here as “lowering computational efficiency” is contradictory; it likely means reducing computational cost or improving efficiency) the computational efficiency of optimization design. Although prior research has enhanced the accuracy and generalization capabilities of RBDO through surrogate modeling, existing studies often involve high-dimensional problems, imposing stricter requirements on surrogate model construction.

Manifold learning methods represent a class of common dimensionality reduction techniques for high-dimensional data spaces.^32,33 Notable algorithms include Laplacian Eigenmaps, ³⁴ locally linear embedding (LLE), ³⁵ diffusion maps (DM), ³⁶ and maximum variance unfolding (MVU), ³⁷ among others. Yuan et al.^38,39 utilized the LE algorithm to project original features into a lower-dimensional space, extracting more representative parameters, and assessed the health status of wind turbines using the standard deviations of horizontal and vertical scales as evaluation metrics. Li and Feng ⁴⁰ combined the LLE algorithm with LE to propose a method that preserves neighborhood information in LLE, employing the “Hotelling’s T² statistic” for industrial process monitoring. Jia et al. ⁴¹ proposed applying DM to equipment degradation assessment modeling and monitoring. In the modeling phase, DM was used to extract eigenvalues for modeling normal equipment behavior, while in the monitoring phase, DM was employed to obtain eigenvalues for each monitoring cycle, enabling degradation assessment through Euclidean distance deviation metrics. Liu et al. ⁴² adopted the MVU algorithm to construct a low-dimensional output space, combined with a Gaussian regression model to capture nonlinear relationships between inputs and outputs, and assessed equipment operational status using the “Hotelling’s T² statistic” and squared prediction error. From existing research, it is evident that projecting a set of high-dimensional data with specific features into a lower-dimensional manifold space using manifold learning can effectively capture the primary information of the original high-dimensional data.

Given the limitations imposed by the substantial time costs associated with constructing surrogate models in the optimization design of complex engineering structures, this study proposes an optimization design strategy for complex engineering structures that incorporates manifold learning. By leveraging manifold learning, the “curse of dimensionality” induced by high-dimensional uncertain variables is mitigated. A data-driven PC-Kriging approach is employed to evaluate limit states for engineering problems, thereby reducing the computational costs incurred by reliability analysis. Finally, a reliability optimization design framework is constructed using a decoupling optimization strategy, abbreviated as L-PCK. The subsequent parts of this study are organized in the following manner: Section 2 introduces the theoretical foundations of the adopted manifold learning strategies and the fundamental principles of optimization design. Section 3 elaborates on the surrogate model construction strategy that integrates manifold learning and provides an overview of the RBDO framework. Section 4 discusses and analyzes the effectiveness and superiority of the proposed method through various case studies. Section 5 concludes the paper and explores potential future applications of the proposed methodology.

Fundamental principles underlying optimization design

Theoretical foundations of manifold learning

In the optimization design process for complex engineering problems, constructing surrogate models based on data is often necessary. However, the presence of numerous uncertainties frequently results in data characterized by large scale, high dimensionality, and strong nonlinearity. To address these challenges, employing dimensionality reduction techniques to simplify data complexity and extract limit state information from complex engineering objects is of significant practical importance. First, for data, that is, structurally complex, high-dimensional, sparse, and noisy, effectively learning its intrinsic structure involves two key aspects. It is essential to capture the local nonlinear characteristics of the data to preserve its complex structure while simultaneously grasping global patterns to provide a crucial basis for decision-making. Second, noise and outliers can adversely affect the stability and accuracy of models during the modeling process. To overcome this difficulty, it is necessary to distinguish the roles of normal and anomalous samples in the modeling process. Additionally, it is imperative to filter core information from the data, eliminate noise and redundancy, achieve an efficient representation of the data, and improve computational efficiency.

Dimensionality reduction of data can be achieved through two primary strategies: feature selection and feature extraction, with the latter also referred to as subspace learning. Feature selection involves selecting the most representative or discriminative features from the original data to form a new feature set, thereby reducing dimensionality while preserving critical information from the original data. Feature extraction, on the other hand, derives a new set of features from the original data using certain criterion functions, or projects the original features into another space. This approach comprehensively considers all inherent information in the data, preventing the loss of useful information. Relevant research has been conducted on reliability analysis and optimization design methods based on dimensionality reduction techniques, with findings demonstrating that this category of methods significantly enhances computational efficiency.^43–45 However, as engineering problems grow increasingly complex, there arises a further demand for optimization design methods rooted in dimensionality reduction techniques. Additionally, the computational cost incurred in establishing surrogate models cannot be overlooked, and reducing the complexity of surrogate model construction strategies is equally of significant importance.

Manifold learning is a type of feature extraction method that excels in dimensionality reduction mapping and modeling. LLE is a representative manifold learning method that assumes data exhibits linearity within local neighborhoods, meaning each data point can be linearly represented by its neighboring points. By preserving this local linear relationship, the algorithm maps high-dimensional data into a low-dimensional space while retaining the original data’s topological structure. Compared to traditional linear dimensionality reduction methods (e.g. PCA), LLE is more effective in handling nonlinear data and uncovering the low-dimensional manifold structures hidden in high-dimensional spaces. Furthermore, LLE achieves dimensionality reduction by solving the eigenvalues of a sparse matrix, resulting in relatively low computational complexity. Below, we introduce a LLE algorithm, which is incorporated into this study.

The core concept of the LLE algorithm ⁴⁶ lies in the assumption that each data point and its neighboring points reside on or near a locally linear patch of the manifold. The local geometric structure of the manifold can be captured by a set of linear reconstruction coefficients, such that the relationships between data samples remain preserved before and after the projection. The main computational procedure of the LLE algorithm is depicted in Figure 1.

Figure 1.

Flowchart of particle swarm optimization algorithm.

Let the high-dimensional dataset $X$ consist of $M$ data points denoted as $X_i$ , where $i=1,2,\dots ,M$ . Each $D_1$ possesses a dimensionality of $X_i$ . The computational procedure can be condensed into the following three steps:

• Step 1: In the LLE algorithm, the $k$ -nearest neighbors algorithm is commonly employed to identify the neighboring points of each sample. However, the neighborhood size $k$ in this algorithm must be predetermined. The choice of $k$ in the LLE algorithm critically determines the extent to which local and global properties of the data are preserved during dimensionality reduction. An excessively large value of $k$ may cause the LLE algorithm to exhibit stronger global characteristics, thereby losing local information between data points. Conversely, an overly small value of $k$ may result in the embedding failing to maintain the spatial topological structure of the original dataset.

• Step 2: Mining the Local Structure of the Original Data Space. Assume that each sample data point $X_i$ is sampled from an underlying manifold. For each data point $X_i$ , we identify its $k$ -nearest neighbors based on Euclidean distance. The weights between local data points are then computed by minimizing the reconstruction error, which is defined as shown in equation (1):

$$\epsilon \left(W\right)=\sum _{i=1}^M\Vert X_i-\sum _{j=1}^k{W}_{\mathrm{ij}}{X}_j{\Vert }_2^2$$

(1)

where $W_{\mathrm{ij}}$ denotes the contribution of the jth local neighbor of the ith data point in reconstructing the ith data point; $X_i$ denotes the jth local neighbor data point of the ith data point. The following two constraints are imposed on $W_{\mathrm{ij}}$ :

• 1. The weights $W_{\mathrm{ij}}$ are non-zero only for the k-nearest neighbors of the ith data point, $j=1,2,\dots k$ .

• 2. The sum of the weights assigned to the ith data point equals 1, as shown in equation (2):

$$\sum _{j=1}^k{W}_{\mathrm{ij}}=1$$

(2)

By optimizing the weights $W_{\mathrm{ij}}$ , retain them as the combination coefficients for computing the low-dimensional embeddings $y_i$ and $y_j$ corresponding to the original data points $x_i$ and $x_j$ . Equation (1) is reformulated into a matrix form as equation (3). The problem is then solved using the Lagrange multiplier method, leading to the establishment of the Lagrangian function as shown in equation (4):

$$\{\begin{array}{c}\epsilon \left(W\right)=\sum _{i=1}^M{W}_i^T{Z}_i{W}_i\hfill \\ \multicolumn{1}{c}{Z_i={(x_i-x_j)}^T(x_i-x_j)}\end{array}$$

(3)

$$L\left(W\right)=\sum _{i=1}^m{W}_i^T{Z}_i{W}_i+\lambda (W_i^T{1}_k-1)$$

(4)

where, $1_k$ stands for a k-dimensional vector with all elements equal to 1. By taking the derivative of equation (4), setting the derivative equal to 0, and normalizing the resulting expression, the following is obtained:

$$W_i=\frac{Z_i^{-1}{1}_k}{1_k^T{Z}_i^{-1}{1}_k}$$

(5)

• Step 3: Obtaining the Low-Dimensional Embedding Results. In the low-dimensional embedding space, we retain the weight matrix computed in Step 2. We then minimize the global reconstruction error function, as shown in equation (6), while incorporating the constraint conditions specified in equation (7). This process yields the standardized low-dimensional embedding results corresponding to the original high-dimensional data:

$$\min \sum _{i=1}^M{\Vert Y_i-\sum _{j=1}^k{w}_{\mathrm{ij}}{Y}_j\Vert }_2^2$$

(6)

$$\{\begin{array}{c}\sum _{i=1}^M{Y}_i=0_d\hfill \\ \multicolumn{1}{c}{\frac{1}{M}\sum _{i=1}^M{Y}_i{Y}_i^T=I}\end{array}$$

(7)

where, $Y_i$ denotes the low-dimensional embedding result corresponding to the original data point $X_i$ ; $Y_j$ represents the jth nearest neighbor of $Y_i$ in the low-dimensional space; and $d$ indicates the dimensionality of the low-dimensional embedded data.

For computational convenience, the objective function can be reformulated as:

$$\begin{array}{c}[Y_1^T{Y}_2^T\dots Y_m^T]=\left[\begin{array}{cccc}{w}_{11}\hfill & w_{12}\hfill & \cdots \hfill & w_{1m}\hfill \\ \multicolumn{1}{c}{w_{21}}& w_{22}\hfill & \cdots \hfill & w_{2m}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{w_{m1}}& w_{m2}\hfill & \cdots \hfill & w_{\mathrm{mm}}\hfill \end{array}\right]\left[\begin{array}{c}{Y}_1\hfill \\ \multicolumn{1}{c}{Y_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{Y_m}\end{array}\right]\hfill \\ =Y{(I-W)}^T(I-W)Y^T=\mathrm{YG}{Y}^T\hfill \end{array}$$

(8)

Based on the relationship between the eigenvalues and eigenvectors of a matrix, we can derive the following:

$$\mathrm{GP}=\mathrm{diag}[{\lambda }_1{\lambda }_2\cdots {\lambda }_d]P$$

(9)

where, $P$ denotes the eigenvector corresponding to the matrix $G$ , and $\lambda$ represents the eigenvalue associated with the matrix $G$ . In general, eigenvectors corresponding to distinct eigenvalues are orthogonal to each other. By multiplying both sides of equation (9) by $P^T$ , obtain:

$$P^T\mathrm{GP}=\mathrm{diag}[{\lambda }_1{\lambda }_2\cdots {\lambda }_d]$$

(10)

From equation (10), it can be observed that the eigenvectors corresponding to the matrix $G$ constitute the solutions to the objective function presented in equation (6). The constraint imposed by equation (7) results in zero being a trivial eigenvalue of the matrix $G$ . In the scenario where the goal is to obtain a $d$ -dimensional low-dimensional embedding, a process is carried out to extract the eigenvectors that correspond to the $d$ smallest non-zero eigenvalues of the matrix $G$ . These extracted eigenvectors then serve as the final embedding results.

Theoretical foundations of optimization design

Evaluation of uncertainty constraints

In the process of RBDO, reliability analysis methods are commonly employed for constraint evaluation. Structural reliability refers to the capability of a structure to maintain its required functions and performance under various external loads, environmental influences, and uncertainties throughout its specified design service life. ⁴⁷ In the course of structural reliability analysis, an evaluation of the limit states is required. Limit states serve as commonly used criteria for assessing whether a structure has failed, and they are governed by a limit-state function (LSF) $G\left(\mathit{x}\right)$ , where $\mathit{x}$ represents random variables:

$$G\left(\mathit{x}\right)=\{\begin{array}{c}>0,\mathrm{safe}\hfill \\ \multicolumn{1}{c}{=0,\mathrm{Limit}\mathrm{state}}\\ \multicolumn{1}{c}{<0,\mathrm{fail}}\end{array}$$

(11)

The joint probability density function of the random variables $\mathit{x}$ is denoted as $f\left(\mathit{x}\right)$ . Thereafter, the failure probability of the structure is determined by:

$$P_f={\int }_{G\left(\mathit{x}\right)\le 0}{f}_p\left(\mathit{x}\right)d\mathit{x}$$

(12)

The reliability is given by:

$$P_r=1-P_f$$

(13)

The reliability $P_r$ is a physical quantity that describes the safety of a structure. A larger value of $P_r$ indicates a safer structure. Conversely, the safety of a structure can also be described by its failure probability $P_r$ , where a smaller value of $P_r$ implies a safer structure. Based on equations (11) and (12), it can be concluded that the failure probability of the structure is equivalent to the probability of the LSF $G\left(\mathit{x}\right)$ being <0. Assuming that the random variable $\mathit{x}$ follows a Gaussian distribution, a standardized transformation is performed on $\mathit{x}$ to convert it into a standard Gaussian variable $\mathit{u}$ , as follows:

$$P_f=P(G\left(\mathit{x}\right)\le 0)=P(g\left(\mathit{u}\right)\le 0)=F_G\left(0\right)$$

(14)

where, $g\left(\mathit{u}\right)$ represents the LSF in the standard normal space, and $F_G$ denotes the cumulative distribution function (CDF) of $g\left(\mathit{u}\right)$ . In engineering practice, the reliability index $\beta$ is also frequently employed to describe the reliability of a structure. The reliability index can be calculated using the following formula:

$$\begin{array}{c}{F}_G\left(0\right)=\mathrm{\Phi }(-\beta )\hfill \\ \beta =-{\mathrm{\Phi }}^{-1}\left(F_G\left(0\right)\right)\hfill \end{array}$$

(15)

where, $\mathrm{\Phi }$ denotes the CDF of a standard normal random variable. Through the application of the Rosenblatt transformation or other appropriate methods, reliability analysis can proceed using the reliability index, despite the random variables exhibiting non-normal distributions and the LSF being nonlinear.

FORM approximates the LSF with a hyperplane tangent to it at the most probable point (MPP). As a result, the reliability index $\beta$ is defined as the distance measured from the origin to the hyperplane that represents the LSF within the standard normal space. Hence, determining the Hasofer–Lind (H-L) reliability index is equivalent to tackling the optimization problem stated below:

$$\begin{array}{c}\min \beta =\Vert \mathit{u}\Vert \hfill \\ \mathrm{s}.\mathrm{t}.g\left(\mathit{u}\right)=0\hfill \end{array}$$

(16)

where, $g\left(\mathit{u}\right)$ represents the limit state surface in the standard normal space, and $\mathit{u}$ is the standardized form of the normally distributed random variable $\mathit{x}$ , defined as follows:

$$\mathit{u}=\frac{\mathit{x}-\mu }{\sigma }$$

(17)

where, $\sigma$ and $\mu$ are the standard deviation and mean of a random variable. The solution ${\mathit{u}}^{*}$ to the optimization equation (16) is designated as the design point or the MPP of failure. For random variables that do not follow a normal distribution, the Rackwitz– Flessler transformation ⁴⁸ must first be applied to convert them into normally distributed random variables. Subsequently, the first two moments of these transformed variables are obtained. Standardization is then carried out according to equation (17), followed by solving the optimization problem equation (16).

Reliability-based optimization design

The RBDO problem can be formulated by introducing probabilistic constraints into the framework of traditional deterministic optimization problems:

$$\begin{array}{c}\mathrm{find}\mathit{d}={[d_1,d_2,\dots ,d_n]}^T\hfill \\ \min f(\mathit{d},\mathit{X})\hfill \\ \mathrm{s}.\mathrm{t}.P_f(G_j(\mathit{d},\mathit{X})\le 0)\le P_{j,t}(j=1,\dots ,N_p)\hfill \\ {\mathit{d}}^L\le \mathit{d}\le {\mathit{d}}^U\hfill \end{array}$$

(18)

where, $\mathit{X}={[X_1,X_2,\dots ,X_m]}^T$ represents the vector of random variables. $G_j(\mathit{d},\mathit{X})$ denotes the jth LSF, and when its value is <0, it indicates that the structure is in a failure state. $P_f(G_j(\mathit{d},\mathit{X})\le 0)$ represents the probability of failure for the jth LSF, while $P_{j,t}$ represents the allowable failure probability for the jth LSF, which is specified arbitrarily. The overall failure probability of the structure can be expressed as:

$$P_f(G_j(\mathit{d},\mathit{X})\le 0)={\int }_{G_j(\mathit{d},\mathit{X})\le 0}{f}_p\left(\mathit{X}\right)\mathrm{dx}=F_{G_j}(0,\mathit{d})$$

(19)

where, $f_p\left(\mathit{X}\right)$ is defined as the joint probability density function (PDF) pertaining to the random variable $\mathit{X}$ , and $F_{G_j}$ represents the CDF of the jth LSF.

The reliability-based optimization problem presented in equation (18) is solvable through the application of the FORM. The RIA and the PMA are two optimization formulations that integrate reliability constraints based on FORM. The computation of failure probability is transformed by them into the determination of the reliability index and the performance measure, respectively. This approach circumvents the need for direct integration of the joint probability density function of random variables. Compared to RIA, PMA exhibits stronger capabilities in handling nonlinear and high-dimensional problems. The core concept of reliability-based optimization design employing the PMA lies in substituting the probabilistic constraint with a constraint on the performance measure of the LSF. The PMA fundamentally represents the inverse of the reliability index method. The optimization formulation equation (18) under the PMA can be expressed as:

$$\begin{array}{c}\mathrm{find}\mathit{d}={[d_1,d_2,\dots ,d_n]}^T\hfill \\ \min f(\mathit{d},\mathit{x})\hfill \\ \mathrm{s}.\mathrm{t}.G_{j,p}(\mathit{d},\mathit{x})\ge 0(j=1,\dots ,N_p)\hfill \\ {\mathit{d}}^L\le \mathit{d}\le {\mathit{d}}^U\hfill \end{array}$$

(20)

where, $G_{j,p}$ represents the performance measure of the LSF $G_j$ , which is calculated using:

$$G_{j,p}(\mathit{d},\mathit{x})=F_{\mathrm{Gj}}^{-1}(\mathrm{\Phi }(-{\beta }_{j,t}),\mathit{d})$$

(21)

The PMA is essentially a constrained optimization problem:

$$\begin{array}{c}\min g\left(\mathit{u}\right)\hfill \\ \mathrm{s}.\mathrm{t}.\Vert \mathit{u}\Vert ={\beta }_t\hfill \end{array}$$

(22)

where, $\mathit{u}$ denotes a standard normal random variable, and ${\beta }_t$ represents the target reliability index. In contrast to the inner optimization problem presented in equation (16) within the PMA, the objective function of the optimization problem defined in equation (22) is the LSF, whereas the constraint function is the 2-norm of the random variables. Thus, the PMA is alternatively termed the inverse reliability analysis method. The inverse reliability analysis method endeavors to find the point which minimizes the LSF’s value and resides on the standard normal spherical surface having a radius equal to the target reliability index. In other words, we aim to find the safest combination of design parameters that meets the reliability index requirement, ensuring that the probability of structural failure or the reliability index attains a preset target value. This point is termed the Minimum performance target point. Figure 2 illustrates the iterative process of the PMA. ${\mathit{u}}_{\mathrm{MPTP}}^{*}$ represents the searched minimum performance target point.

Figure 2.

Schematic diagram of the iterative process of the PMA.

In this study, the first-order inverse reliability analysis method adopts the hybrid mean value method proposed by Youn et al. ⁴⁹ This method is capable of distinguishing the convexity or concavity of the limit state function and subsequently solving for the Most Probable Target Point, making it particularly suitable for nonlinear and high-dimensional problems. However, since the method involves direct evaluation and optimization of performance measures, it may consume substantial computational resources and time. To address this, the present study enhances the computational efficiency of the method by incorporating surrogate modeling techniques.

Surrogate model-assisted RBDO

Data-driven surrogate model construction

During the process of optimization design, it is often necessary to repeatedly evaluate the limit states of constraints multiple times. To enhance the computational efficiency and accuracy for complex models, various surrogate models have been developed to replace time-consuming physical simulations or computations involving implicit functions. PCE and the Kriging model are commonly used surrogate models. In this study, we employ a combined surrogate model derived from these two models, which demonstrates higher computational accuracy and efficiency compared to using either model individually.

The Kriging model is classified as a stochastic algorithm. It assumes the model output $M^K\left(\mathit{x}\right)$ to be a realization of a Gaussian process, with the index $\mathit{x}\in \mathrm{DX}\subset \mathrm{RM}$ , and is described by the equation provided below:

$$M^K\left(\mathit{x}\right)={\mathit{b}}^T\mathit{f}\left(\mathit{x}\right)+{\sigma }^{2}Z(\mathit{x},\omega )$$

(23)

where, the first term ${\mathit{b}}^{T}f\left(\mathit{x}\right)$ denotes the Gaussian process’s mean, which is composed of P arbitrary functions $\{f_j;j=1,\dots ,P\}$ and corresponding coefficients $\{{\beta }_j;j=1,\dots ,P\}$ . The second term consists of the variance of the Gaussian process, denoted as ${\sigma }^2$ , and a stationary Gaussian process $Z(\mathit{x},\omega )$ with zero mean and unit variance. The fundamental probability space is symbolized by $\omega$ and is characterized by the correlation function $R$ along with its hyperparameter $\theta$ . The correlation function $R=R(\mathit{x},\mathit{x}\prime ;\mathbf{\theta })$ describes the correlation between two sample points in the output space and depends on $\mathit{x}$ , $\mathit{x}\prime$ , and the hyperparameters $\mathbf{\theta }$ .

PCE ⁵⁰ serve as a surrogate modeling approach that leverages polynomial chaos to approximate the stochastic responses of complex systems. The fundamental idea behind PCE is to utilize the properties of orthogonal polynomials to approximate the probability distribution of the stochastic response. Denote the set of orthogonal basis functions as $\{\psi \alpha :\alpha \in N^d\}$ . Consider $\mathit{x}$ to be a d-dimensional random vector with independent components, defined on a domain $D\in R^d$ , where the probability density function along each dimension is square-integrable. Then, the following holds:

$$y=M\left(x\right)\approx M^{\mathrm{PCE}}\left(\mathit{x}\right)=\sum _{\alpha \in A}{c}_{\alpha }{\psi }_{\alpha }\left(\mathit{x}\right)$$

(24)

where, $M$ represents the original model, $M^{\mathrm{PCE}}$ denotes the surrogate model constructed using PCE, $y$ signifies the stochastic response, and $A\in N^d$ stands for a finite set of orthogonal basis functions. The coefficients $c_{\alpha }$ can be determined through regression methods. Let ${\left\{{\mathit{x}}^{\left(k\right)}\right\}}_{k=1}^N\subset D$ be a sample set in the input space, and $\mathit{Y}={(y^{\left(1\right)},y^{\left(2\right)},\dots ,y^{\left(N\right)})}^T$ be a vector comprising the model responses $y^{\left(k\right)}=M\left({\mathit{x}}^{\left(k\right)}\right)$ . The regression matrix $\mathit{\Psi }$ is defined based on ${\Psi }_{\mathrm{ij}}={\psi }_j\left(x^{\left(i\right)}\right)$ . By solving the following optimization problem:

$$\hat{\mathit{c}}=\mathrm{aig}\underset{c\in R^d}{\min }\Vert {\mathit{\Psi }}_c-\mathit{Y}{\Vert }_2$$

(25)

Owing to the computational cost determined by the total number of terms in the expansion, the coefficient polynomial expansion method is utilized for solving high-dimensional problems. By introducing the L1-norm, the problem of solving for the coefficients can be formulated as follows:

$$\hat{\mathit{c}}=\arg \underset{\mathit{c}}{\min }\Vert {\mathit{\Psi }}_{\mathit{c}}-\mathit{Y}{\Vert }_2^2+\lambda \Vert \mathit{c}{\Vert }_1$$

(26)

where, the L1-norm is represented with $\Vert \mathit{c}{\Vert }_1=\sum _i{c}_i$ , and $\lambda$ serves as the Lagrange multiplier. The approximate solution to the L1-norm minimization problem is obtained using the Least Angle Regression algorithm.

To further enhance the accuracy and efficiency of surrogate models in substituting for physical simulations, the PC-Kriging model is put forward, drawing on the principles of Kriging and PCE. Essentially, it is an improved Kriging model that employs sparse polynomial chaos bases as the regression basis functions. The advantage of this model lies in its ability to capture the global behavior of the computational model using regression-type PCE while also capturing local variations through interpolation-type Kriging. The basis functions of this model are described as follows:

$$y\approx M^{\mathrm{PCK}}\left(\mathit{x}\right)=\sum _{\alpha \in A}{y}_{\alpha }{\psi }_{\alpha }\left(\mathit{x}\right)+{\sigma }^{2}Z(\mathit{x},\omega )$$

(27)

where, $\sum _{\alpha \in A}{y}_{\alpha }{\psi }_{\alpha }\left(\mathit{x}\right)$ represents the basis function of the PC-Kriging model, which is expressed as a weighted sum of standard orthogonal polynomials. Similar to the Kriging model, ${\sigma }^2$ and $Z(\mathit{x},\omega {)}^2$ denote, respectively, a variance parameter along with a stationary Gaussian random process featuring zero mean and unit variance is considered. The predicted value is obtained from equation (27).

Reliability-based design optimization using surrogate models

A new RBDO framework can be established by substituting the limit state evaluation in the optimization problem with the constructed surrogate model. To further reduce computational costs, the construction of an adaptive PC-Kriging model is incorporated. The procedures for computation are outlined below:

• Step 1: The initial candidate basis for the basis functions is a set of orthogonal polynomials with orders $p=p_0$ and q-norm $q=0.5$ . The inputs $\mathit{X}$ for the initial sample points are obtained through Latin Hypercube Sampling (LHS), and their corresponding response values $\mathit{Y}$ are computed. Let the number of initial sample points be $N=N_0$ , where:

$$\begin{array}{c}{\mathit{\Omega }}_0=\{(x_i,y_i),i=1,\dots ,N_0\},\hfill \\ \mathit{X}=\{x_1,x_2,\dots ,x_{N_0}\}\hfill \end{array}$$

(28)

• Step 2: Calculate the Leave-One-Out (LOO) error of the model:

$${\epsilon }_{\mathrm{LOO}}=\frac{\sum _{i=1}^N{(M\left(x_i\right)-M^{\mathrm{PCK}\backslash i}\left(x_i\right))}^2}{\sum _{i=1}^N{(M\left(x_i\right)-{\hat{\mu }}_{\mathit{Y}})}^2}$$

(29)

where, $M$ represents the surrogate model formed by the existing set of sample points, $M^{\mathrm{PCK}\backslash i}$ is the model generated from the sample set after excluding the ith sample point, and ${\hat{\mu }}_Y$ is the mean of the sample point responses $Y$ . The convergence criterion is based on whether the LOO error of model exceeds a predefined maximum error threshold ${\epsilon }_{\max }$ . If the LOO error of $M$ is higher than ${\epsilon }_{\max }$ , the convergence criterion is not met, and the process proceeds to Step 3. If the criterion is satisfied, the iterative process is terminated, yielding the final surrogate model $M$ .

• Step 3: Modify the candidate basis for the basis functions. If $p=1$ , then set $p=p+1$ . If $p>1$ , first increase the q-norm by setting $q=q+0.25$ . When $q$ equals 1, increment the polynomial order by setting $p=p+1$ and reset $q$ to 0.5.

• Step 4: Draw $N_{\mathrm{add}}$ new sample points using the LHS method and compute their corresponding responses. Add these points to the existing set of sample points. The new set of sample points becomes:

$$\begin{array}{c}{\mathit{\Omega }}_{\mathrm{add}}=\{(x_j,y_j),j=1,\dots ,N_{\mathrm{add}}\},\hfill \\ \mathit{\Omega }={\mathit{\Omega }}_0\cup {\mathit{\Omega }}_{\mathrm{add}},\hfill \\ \mathit{X}=\{x_0,x_1,\dots ,x_{N_0},x_{N_0+1},\dots ,x_{N_0+N_{\mathrm{add}}}\}\hfill \end{array}$$

(30)

A new PC-Kriging model $M$ is generated based on the updated set of sample points. The LOO error of the model is calculated according to equation (29). If the convergence criterion is not met, the process proceeds to Step 3. If the criterion is satisfied, the iterative process is terminated, and the final surrogate model $M$ is obtained.

The proposed adaptive PC-Kriging model is employed to replace the original model for reliability analysis and optimization design. The computational procedure for RBDO using the surrogate model is illustrated in Figure 3.

Figure 3.

Reliability-based optimization process based on the PC-Kriging model.

Optimization design framework incorporating manifold learning

An improved method for constructing surrogate models

As described in Section 2.3.2, while it is feasible to obtain the desired PC-Kriging model based on data-driven approaches, complex engineering problems often encounter challenges due to the high dimensionality of parameters. To address this issue, this study proposes an improved method for constructing surrogate models by incorporating dimensionality reduction techniques. The steps for constructing the new surrogate model are described as follows:

• Step 1: Clearly define the optimization objective and the constraint functions. Construct a high-dimensional dataset $\left\{{\mathit{X}}_m\right\}$ based on the black-box problems within the constraint functions. The total number of sample points in the dataset is m, and assume the original dimensionality of the data is $n_{\mathrm{high}}$ .

• Step 2: Apply the LLE strategy to perform dimensionality reduction on the constructed high-dimensional dataset, thereby obtaining a low-dimensional dataset $\left\{{\mathit{X}}_{m1}\right\}$ . Assume the reduced dimensionality is $n_{\mathrm{low}}$ .

• Step 3: Train the PC-Kriging model based on the training dataset $\left\{{\mathit{X}}_{m1}\right\}$ that has undergone dimensionality reduction.

• Step 4: Calculate the LOO error $\epsilon$ of the PC-Kriging model.

• Step 5: Determine whether the error meets the stopping criterion ${\epsilon }_{\max }=1\times {10}^{-3}$ . If it does, output the current surrogate model; otherwise, proceed to Step 6.

• Step 6: Regenerate a candidate sample set $\left\{{\mathit{X}}_{\mathrm{lhs}}\right\}$ using the LHS method, and apply the LLE strategy to perform dimensionality reduction on the entire candidate sample set $\left\{{\mathit{X}}_{\mathrm{lhs}}\right\}$ . Based on the theory outlined in Section 2.3.2, resample new sample points, and add these new sample points $\left\{{\mathit{X}}_{\mathrm{madd}}\right\}$ to the collected low-dimensional dataset $\left\{{\mathit{X}}_{m1}\right\}$ . The proposed strategy for constructing the surrogate model is illustrated in Figure 4.

Figure 4.

Flowchart for the construction of an adaptive PC-Kriging model assisted by manifold learning.

Overall optimization design framework

By coupling the surrogate model construction process with optimization design, the RBDO framework proposed in this paper is ultimately obtained. During the optimization process, a heuristic optimization algorithm is employed to address the problem presented in equation (20). The specific workflow of the new optimization framework is as follows:

• Step 1: Define the optimization objectives and constraint functions, and establish the optimization model.

• Step 2: Generate an initial candidate sample set and construct the initial PC-Kriging model based on the methodology described in Section 3.1.

• Step 3: Replace the time-consuming constraint function evaluations with the constructed PC-Kriging model, and reconstruct the surrogate model-based optimization model.

• Step 4: Based on the constructed optimization model, transform the probabilistic constraints into deterministic constraints using the PMA as described in equation (22).

• Step 5: Establish the current fitness function based on the deterministic constraint optimization problem and compute the current optimal design point.

• Step 6: Check whether the changes in the objective function and constraints meet the specified tolerance. If the constraints are satisfied, output the optimal design point and the objective function value; otherwise, return to Step 4 for iterative optimization. Figure 5 presents the flowchart of the overall RBDO process.

Figure 5.

Overall optimization design framework based on manifold learning.

The deterministic optimization algorithm employed in this study is a modified Particle Swarm Optimization algorithm enhanced with a Harmony Search-based updating strategy. ⁵¹ The computational procedure of this algorithm is outlined as follows:

Step number	Execute operation
1.	Parameter initialization: number of particles, maximum iteration count, inertia weight, acceleration coefficients, harmony search coefficient, dynamic adjustment rate
2.	Initialize particle swarm and harmony memory
3.	Update inertia weight and particle states
4.	Hybrid PSO-harmony search position update with boundary constraints
5.	Convergence check. If converged: Output global best solution; else: return to Step 3

Case study analysis

Mathematical examples

A 10-dimensional nonlinear optimization numerical example

The optimization model for this case study is presented as follows ⁵² :

$$\{\begin{array}{c}\min f\left(\mathit{d}\right)=d_1^2+d_2^2+d_1{d}_2-14d_1-16d_2+\hfill \\ {(d_3-10)}^2+4{(d_4-5)}^2+{(d_5-3)}^2+\hfill \\ 2{(d_6-1)}^2+5d_7^2+7{(d_8-11)}^2+\hfill \\ 2{(d_9-10)}^2+{(d_{10}-7)}^2+45\hfill \\ \mathrm{s}.\mathrm{t}.P_f(G_i\ge 0)\ge \mathrm{\Phi }(3.0)i=1,2,\dots ,8\hfill \\ 0\le d_i\le 15i=1,2,\dots ,10\hfill \end{array}$$

(31)

where, $\mathit{d}$ represents the mean value of the random variable $x$ . $P_f(·)$ denotes the probability of event occurrence. $G_i(·)$ represents distinct limit state functions. The LSF for this example is expressed as:

$$G_1\left(\mathit{x}\right)=1-\frac{4x_1+5x_2-3x_7+9x_8}{105}$$

(32)

$$G_2\left(\mathit{x}\right)=-10x_1+8x_2+17x_7-2x_8$$

(33)

$$G_3\left(\mathit{x}\right)=1+\frac{8x_1-2x_2-5x_9+2x_{10}}{12}$$

(34)

$$G_4\left(\mathit{x}\right)=1-\frac{3{(x_1-2)}^2+4{(x_2-3)}^2+2x_3^2-7x_4}{120}$$

(35)

$$G_5\left(\mathit{x}\right)=1-\frac{5x_1^2+8x_2+{(x_3-6)}^2-2x_4}{40}$$

(36)

$$G_6\left(\mathit{x}\right)=1-\frac{0.5{(x_1-8)}^2+2{(x_2-4)}^2+3x_5^2-x_6}{120}$$

(37)

$$G_7\left(x\right)=-x_1^2-2{(x_2-2)}^2+2x_1{x}_2-14x_5+6x_6$$

(38)

$$G_8\left(\mathit{x}\right)=3x_1^2-6x_2-12{(x_9-8)}^2+7x_{10}$$

(39)

where, the optimization variables $x_1$ to $x_{10}$ all follow a normal distribution, with each optimization variable having a mean value of $\mu ={\{2.17,2.36,8.77,5.10,0.99,1.43,1.32,9.83,8.28,8.38\}}^T$ and a standard deviation of 0.02.

As shown in Table 1, this study compares the proposed method with deterministic optimization approaches, uncertainty optimization methods based on single surrogate models (including polynomial chaos expansion (PCE), Kriging, and sparse polynomial chaos expansion (SPCE)), and uncertainty optimization methods employing alternative dimensionality reduction strategies (e.g. HDRA ⁴³ ).

Table 1.

Comparison of optimization results from different methods.

Method	Number of sample points	Optimal design point	Optimal objective function	LOO error	Iteration count
No model	–	2.1385, 2.3298, 8.7092, 5.1018, 0.9232, 1.4589, 1.3869, 9.8071, 8.1564, 8.4618	27.7226	–	6
PCE	2079	2.1354, 2.3219, 8.7109, 5.0972, 0.9295, 1.4617, 1.3892, 9.8083, 8.1595, 8.4612	27.7766	2.899 × 10−3	10
SPCE	1566	2.1327, 2.3315, 8.7094, 5.1010, 0.9293, 1.4642, 1.3890, 9.8097, 8.1574, 8.4681	27.7301	3.6344 × 10−3	8
Kriging	1683	2.1348, 2.3316, 8.7086, 5.1024, 0.9247, 1.4527, 1.3879, 9.8077, 8.1567, 8.4768	27.7637	4.2303 × 10−3	8
HDRA	544	2.1358, 2.3217, 8.7082, 5.1010, 0.9291, 1.4624, 1.3850, 9.8097, 8.1527, 8.4658	27.7705	1.6540 × 10−3	5
L-PCK	435	2.1368, 2.3319, 8.7082, 5.1013, 0.9277, 1.4603, 1.3899, 9.8101, 8.1532, 8.4632	27.7213	1.7356 × 10−3	5

As can be seen from the aforementioned results, the optimization result obtained without substituting the constraint function with a surrogate model is 27.7226. By comparing the computational results assisted by different surrogate models, it is evident that the optimization result derived from L-PCK exhibits the smallest error. Additionally, when comparing the final LOO errors of different surrogate models, the L-PCK model demonstrates higher accuracy, corroborating the precision of the optimization result. The LOO errors of L-PCK and HDRA exhibit no significant difference. Moreover, when comparing the total number of sample points utilized by various optimization methods, it becomes evident that L-PCK entails the lowest computational cost, thereby highlighting the superiority of the proposed approach.

For the purpose of visually demonstrating the variations in the objective function value across each iteration, the iterative curve of the objective function corresponding to this case study is presented in Figure 6.

Figure 6.

The iterative history of the objective function.

An 11-dimensional vehicle crash case study

The vehicle side-impact problem is one of the widely adopted standard case studies in the field of RBDO. ⁵³ The initial optimization model consists of nine design random variables, two non-design random variables, and 10 constraint conditions However, only three of these constraints are effective. In this study, we consider precisely the same three effective constraints as those in Jung et al. ⁵³ :

$$\begin{array}{c}\min f\left({\mathit{\mu }}_x\right)=1.98+4.9d_1+6.67d_2+6.98d_3\hfill \\ +4.01d_4+1.78d_5+2.73d_7\hfill \\ \mathrm{s}.\mathrm{t}.P_f(G_j\left(x\right)\ge 0)\ge \mathrm{\Phi }(1.28),j=1,2,3\hfill \\ {\mathit{d}}^L\le \mathit{d}\le {\mathit{d}}^U,x_i~N({\mu }_i,{\sigma }_i^2),i=1,2,\dots ,11\hfill \\ \mathrm{where}\hfill \\ G_1=14.35-9.9x_2-12.9x_1{x}_8+0.1107x_3{x}_{10}\hfill \\ G_2=0.72-0.5x_4-0.19x_2{x}_3-0.0122x_4{x}_{10}\hfill \\ +0.009325x_6{x}_{10}+0.000191x_{11}^2\hfill \\ G_3=1.35-0.489x_3{x}_7-0.843x_5{x}_6+0.0432x_9{x}_{10}\hfill \\ -0.0556x_9{x}_{11}-0.000786x_{11}^2\hfill \\ {\mathit{d}}^L=[0.5,0.45,0.5,0.5,0.875,0.4,0.4,0.192,0.192]\hfill \\ {\mathit{d}}^U=[1.5,1.35,1.5,1.5,2.625,1.2,1.2,0.345,0.345]\hfill \end{array}$$

(40)

where, $x_1~x_9$ represent design random variables that follow independent normal distributions, while $x_{10}$ and $x_{11}$ are non-design random variables. Although $x_{10}$ and $x_{11}$ exhibit randomness, their mean values remain unchanged during the optimization process, meaning they are not subject to design adjustments. The sample means and sample standard deviations are given as $[1,1,1,1,2,1,0.3,0.3,0.3,1,-1]$ and $[0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.006,0.06,10,10]$ , respectively.

The optimization results obtained through different approaches are displayed in Table 2.

Table 2.

Comparison of optimization results from different methods.

Method	Number of sample points	Optimal design point	Optimal objective function	LOO error	Iteration count
No model	–	0.5, 1.1442, 0.7684, 1.5, 1.4697, 1.2, 0.4, 0.3127, 0.345	27.1483	–	13
PCE	189	0.5, 1.1189, 0.8729, 1.5, 1.6092, 1.2, 0.4, 0.3095, 0.345	27.9573	1.5311 × 10−4	18
SPCE	212	0.5, 1.1127, 0.8099, 1.5, 1.6376, 1.2, 0.4, 0.3363, 0.345	27.5267	5.0911 × 10−3	16
Kriging	126	0.5, 1.1376, 0.7998, 1.5, 1.5367, 1.2, 0.4, 0.3226, 0.345	27.4427	1.2357 × 10−4	15
HDRA	100	0.5099, 1.1132, 0.8592, 1.5, 1.6351, 1.2, 0.4, 0.3368, 0.345	27.9182	1.1147 × 10−5	13
L-PCK	89	0.5133, 1.0213, 0.8637, 1.5, 1.6234, 1.2, 0.4, 0.3448, 0.345	27.3322	1.3462 × 10−5	12

As can be observed from the table, L-PCK demonstrates higher efficiency and accuracy. A comparison between the optimization design results based on a single surrogate model and those based on the proposed optimization algorithm reveals that L-PCK yields a smaller objective function value. Compared with HDRA implementations that rely on alternative dimension-reduction strategies, L-PCK demonstrates superior computational-resource efficiency in the present numerical example. Moreover, upon comparing the errors of various surrogate models, it becomes clear that the model derived via L-PCK displays the least error, thus emphasizing the superiority of the proposed approach.

In order to more comprehensively illustrate the changes in the objective function throughout the iterative process, Figure 7 displays the numerical values of the objective function derived from the iterative calculations of the proposed method. As can be observed from the figure, L-PCK achieved stable and accurate results after 12 iterations.

Figure 7.

The iterative history of the objective function.

A framework structure case study

This study investigates the reliability analysis and optimization design of a framework structure, ⁵⁴ as illustrated in Figure 8. The numbers enclosed in circles and parentheses correspond to node numbers and element numbers, respectively. All beams in the structure are modeled as Euler-Bernoulli beams with hollow square cross-sections, where the outer side length is denoted as $D_i,i=1,2,3$ and the thickness as $0.05D_i$ . For beam elements $(1,2,3),(4,5,6,7,8,9),(10,11,12)$ , the outer side lengths are specified as $D_1$ , $D_2$ , and $D_3$ , respectively. The constraints are imposed on the horizontal and vertical displacements of Node 1, as well as the vertical displacement of Node 4. The loading conditions consist of concentrated loads $P_1$ and $P_2$ acting in the y-direction at Nodes 2 and 3, along with linearly distributed loads defined by $q_1,\dots ,q_6$ acting on Elements 1, 2, and 3. The structural uncertainty is characterized by 14 parameters $\mathbf{\xi }=({\xi }_1,{\xi }_2,\dots ,{\xi }_{14})$ , each following a uniform distribution over the interval $[-1,1]$ . These parameters correspond to the random expressions of the elastic modulus $E$ , cross-sectional side length $\mathit{D}$ , nodal coordinates ${\mathit{N}}_{\mathrm{xy}}$ , and loads $\mathit{F}$ as follows:

$$E=1000(1+0.2{\xi }_1)$$

(41)

$$\begin{array}{c}\mathit{D}=(D_1,D_2,D_3)\hfill \\ =(d_1+0.05{\xi }_2,d_2+0.05{\xi }_3,d_3+0.05{\xi }_4)\hfill \end{array}$$

(42)

$${\mathit{N}}_{\mathrm{xy}}=\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{1+0.0910{\xi }_5}& 0.0963{\xi }_6\hfill \\ \multicolumn{1}{c}{2+0.0910{\xi }_5}& 0.0971{\xi }_6\hfill \\ \multicolumn{1}{c}{3+0.0699{\xi }_5}& 0\hfill \\ \multicolumn{1}{c}{0.5+0.0868{\xi }_5}& 0.866+0.0684{\xi }_6\hfill \\ \multicolumn{1}{c}{1.5+0.0938{\xi }_5}& 0.866+0.0832{\xi }_6\hfill \\ \multicolumn{1}{c}{2.5+0.0868{\xi }_5}& 0.866+0.0985{\xi }_6\hfill \end{array}\right]$$

(43)

$$\begin{array}{c}\mathit{F}=(P_1\frac{{\xi }_7+1}{2},P_2\frac{{\xi }_8+1}{2},q_1\frac{{\xi }_9+1}{2},\dots ,q_6\frac{{\xi }_{14}+1}{2})\hfill \\ =(25+25{\xi }_7,25+25{\xi }_8,50+50{\xi }_9,\dots ,50+50{\xi }_{14})\hfill \end{array}$$

(44)

Figure 8.

Schematic diagram of a framework structure.

To validate the accuracy and efficiency of the proposed PC-Kriging model, functional relationships between stochastic responses and random variables were constructed using different surrogate models. Where, $d_i=0.6,i=1,2,3$ was adopted. The stochastic responses considered were the vertical displacement of Node 2 and the overall compliance of the truss structure. For the PCE and SPCE models, the polynomial order was set to 4. Within the proposed method, the initial count of sample points, denoted as $N_0$ ., was established at 100, and an extra 50 sample points, labeled as $N_{\mathrm{add}}$ , were incorporated at each iteration step. The maximum LOO error ${\epsilon }_{\max }$ was set to $5\times {10}^{-4}$ . LHS was employed for all sampling processes. The reliability analysis results based on Monte Carlo simulation (MCS) with ${10}^6$ samples from the original model were used as the benchmark. Comparisons of reliability analysis results for displacement response and structural compliance across different surrogate models are presented in Tables 3 and 4, respectively. The iterative history of fitting errors during the construction of the proposed surrogate model is illustrated in Figure 9.

Table 3.

Comparison of reliability analysis results for displacement response using different surrogate models.

Methods	Samples point	Number of polynomial terms	Mean value	Standard deviation	LOO fitting error
MCS	1 × 106	–	3.3492	0.9083	–
PCE	3000	3040	3.2909	0.9727	9.7571 × 10−8
SPCE	500	243	3.3480	0.9035	1.8891 × 10−4
Kriging	500	–	3.3472	0.8972	1.5264 × 10−3
L-PCK	300	149	3.3489	0.9080	6.7810 × 10−5

Table 4.

Comparison of reliability analysis results of different surrogate models for compliance response.

Methods	Samples point	Number of polynomial terms	Mean value	Standard deviation	LOO fitting error
MCS	1 × 106	–	561.6326	250.8855	–
PCE	3000	3040	561.8516	250.0960	3.4952 × 10−7
SPCE	500	243	560.9066	249.1014	2.2027 × 10−4
Kriging	500	–	560.6617	247.7832	6.3985 × 10−4
L-PCK	300	149	561.9304	250.1728	3.7308 × 10−4

Figure 9.

(a) The Loo error of the LPCK model for fitting displacement response during the iteration process; (b) The Loo error of the L-PCK model for fitting compliance response during the iteration process.

As can be seen from tables and Figure 9, the proposed surrogate model construction strategy offers higher accuracy while incurring lower computational costs.

Reliability-based optimization (RBO) design is considered for the truss structure, with the objective of minimizing the total structural volume while satisfying a maximum compliance constraint. The optimization formulation is expressed as follows:

$$\begin{array}{c}\underset{\mathit{d}}{\min }V\hfill \\ \mathrm{s}.\mathrm{t}.P_f[C\ge \overline{C}]\ge \mathrm{\Phi }(1.28)\hfill \\ (0.5,0.5,0.5)\le d\le (1.5,1.5,1.5)\hfill \end{array}$$

(45)

Where, $C$ denotes the computed maximum flexibility response. The maximum compliance $\overline{C}$ is set to 500. When the flexibility response exceeds the predefined maximum allowable flexibility, the structure is considered to have failed. The initial values of the design variables are set to $[0.6,0.6,0.6]$ , and the corresponding deterministic structural optimization formulation is expressed as follows:

$$\begin{array}{c}\underset{\mathit{d}}{\min }V\hfill \\ \mathrm{s}.\mathrm{t}.C-\overline{C}<0\hfill \\ (0.5,0.5,0.5)\le \mathit{d}\le (1.5,1.5,1.5)\hfill \end{array}$$

(46)

The RBDO is performed by replacing the truss structure’s compliance function with the proposed adaptive PC-Kriging model and other surrogate models, respectively. The initial number of sample points $N_0$ is set to 100, with an additional 50 sample points $N_{\mathrm{add}}$ added at each iteration step. The maximum LOO error ${\epsilon }_{\max }$ is set to ${10}^{-3}$ . The polynomial order for the SPCE model is set to 4. A comparison of the optimization results is presented in Table 5.

Table 5.

Comparison of optimization results from different methods.

Method	Number of sample points	Optimal design point	Optimal objective function	LOO error	Iteration count
DO	–	0.7390, 0.5, 0.5	0.6915	–	7
PCE	237	0.9983, 0.5, 0.5	0.9488	1.6231 × 10−2	11
SPCE	198	1.001, 0.5, 0.5	0.9512	4.9888 × 10−4	8
Kriging	253	1.0009, 0.5, 0.5	0.9510	4.4404 × 10−3	8
HDRA	178	1.0009, 0.5, 0.5	0.9510	1.5499 × 10−3	8
L-PCK	150	1.0009, 0.5, 0.5	0.9510	1.6978 × 10−3	8

As can be observed from the results in the table, the method proposed in this paper consumes the least computational cost. Additionally, the L-PCK model exhibits the smallest error, yielding accurate optimization results.

To further illustrate the evolution of the objective function during the iterative process, Figure 10 presents the numerical values of the objective function obtained using the proposed method at each iteration step. As can be observed from the figure, the proposed method demonstrates notable advantages over alternative approaches in terms of both efficiency and accuracy, while exhibiting excellent convergence behavior.

Figure 10.

Iteration of the optimization design of the truss structure.

Engineering case study

Offshore wind turbine support structures

This case study addresses a RBO design problem for a jacket-type offshore wind turbine. ⁵⁵ Figure 11 depicts a geometric schematic of the offshore wind turbine support structure, specifically tailored for a 5 MW wind turbine situated in near- to mid-offshore areas. As illustrated in the figure, the support structure comprises four legs connected to a grouted foundation at the seabed.

Figure 11.

(a) Presents the schematic diagram of the engineering actual conditions of the offshore wind turbine; (b) Presents the schematic diagram of the offshore wind turbine support structure with geometric properties.

The geometric properties of each jacket member are described in Table 6, where r denotes the diameter of the tube, and $\mathit{t}$ signifies the thickness of the tube. The geometric interpretations of the variables, as indicated by their subscripts, are illustrated in Figure 11(b).

Table 6.

Property parameters in this example.

Variables	Min (m)	Max (m)
D_top	0.8	1.3
t_top	0.025	0.045
t_leg1	0.020	0.040
t_leg2	0.025	0.045
t_leg3	0.030	0.050
t_leg4	0.035	0.055
t_leg5	0.045	0.065
D_chord	0.3	0.9
t_chord1	0.008	0.020
t_chord2	0.010	0.025

The finite element model of the offshore wind turbine employs PIPE59 elements to model the beam structures and SOLID95 elements to construct the bottom connection transitions. The materials used in the model are high-strength steel and concrete, with their respective properties listed in Table 7. Figure 12 illustrates the meshing results obtained during the finite element analysis (FEA) of the offshore wind turbine.

Table 7.

Material properties in this example.

Material	Steel	Concrete
Elastic modulus, $E\left(\mathrm{GPa}\right)$	210	70
Poisson’s ratio, $\nu$	0.30	0.19
Density, $\rho (\mathrm{kg}/{\mathrm{m}}^3)$	7850	2000
Buckling strength, $f_y\left(\mathrm{MPa}\right)$	355	-
Compressive strength, $f_c\left(\mathrm{MPa}\right)$	-	200
Tensile strength, $f_t\left(\mathrm{MPa}\right)$	-	10

Figure 12.

Meshing results.

The operational conditions of offshore wind turbines are complex, as they are subjected to various loading conditions that must be considered during the design process.^56,57 In the validation of this case study, multiple uncertain loads were accounted for, among which the dead loads imposed by the superstructure on the support structure are presented in Table 8.

Table 8.

Material properties in this example.

Component	Mass
Transition piece	666,000 kg
Tower	347,460 kg
Rotor nacelle assembly	350,000 kg
Blades	17,740 kg

In addition, under extreme environmental conditions, wind forces can exert a significant impact on offshore wind turbines. ⁵⁸ The loads induced by wind include aerodynamic loads, which are composed of thrust forces, overturning moments, and torsional moments. To account for these effects, this case study also considers the influence of aerodynamic loads by simulating the worst-case scenario where all loads are applied to the jacket side in the global X-direction. The aerodynamic loads of the wind turbine include a thrust force of 781 kN, a tilting moment of 38,567 kN × m, and a torsional moment of 7876 kN × m.

The magnitude of wind loads imposed on this support structure is contingent upon the wind speed. The wind speed measured at the nacelle is typically used as a reference for determining the turbine’s operational speed, while wind speeds at other elevations are calculated using the following equation ⁴² :

$${\overline{V}}_h={\overline{V}}_r{\left(\frac{h}{h_r}\right)}^a$$

(47)

where $\overline{V}$ denotes the velocity of different height. ${\overline{V}}_r$ denotes the wind speed of cabin height. $h$ denotes the height of wind turbine. $h_r$ denotes the height of wind cabin. $a=0.115$ , denotes the roughness factor.

The wind pressure $P_{\mathrm{wind}}$ on the support structure depends on the wind speed ${\overline{V}}_h$ , as expressed by the following equation ⁵⁵ :

$$P_{\mathrm{wind}}=0.613\times {10}^{-3}{\overline{V}}_h^2$$

(48)

The magnitude of the current load under steady flow conditions is given by the following equation:

$$F_{\mathrm{current}}=C_D\frac{\rho }{2}{v}_c^{2}A$$

(49)

where $F_{\mathrm{current}}$ denotes the current load, $\rho$ denotes the density of sea water, $C_D$ denotes the drag coefficient, $A$ denotes the area of model in a plane perpendicular to the flow direction, $v_c$ denotes the current velocity.

Here’s how to calculate wave loads $F_{\mathrm{wave}}$ using Morison’s equation:

$$F_{\mathrm{wave}}=\frac{1}{2}\rho C_D\left|w_n\right|w_D+\frac{1}{4}{\rho }_w{C}_M\pi D^{2}w{\prime }_n$$

(50)

where $w_n$ denotes the horizontal velocity of water particles perpendicular to the axis of the component, $w{\prime }_n$ is the first derivative of $w_n$ . $C_M$ is the inertia force coefficient, taken as 1.0 and 1.6. The delineations of loads and constraints for the model of the offshore wind turbine support structure are presented in Figure 11(a).

Best design for offshore wind turbine support structures

Due to the manufacturing costs of offshore wind turbine support structures and to make them stiffer, the math model for this example is as in equation (51). In this model, the objective function for the optimization design is set as the total mass of steel required, while the constraint functions are defined as stress constraints, deformation constraints, and stability requirements:

$$\begin{array}{c}\min f(\mathit{D},\mathit{t})=M(\mathit{D},\mathit{t})\hfill \\ \mathrm{s}.\mathrm{t}.\mathrm{Pr}(g_i(\mathit{D},\mathit{t})\ge 0)\ge \mathrm{\Phi }(3.719),i=1,...,3\hfill \\ g_1=\left[\sigma \right]-\sigma \hfill \\ g_2=\left[u\right]-u_{max}\hfill \\ g_3=\frac{50E}{f_y}-\frac{\mathit{D}}{\mathit{t}}\hfill \end{array}$$

(51)

where $\mathit{D}$ and $\mathit{t}$ denote the design variables, their value ranges are shown in the Table 6. $M(·)$ denotes the structure volume. $\left[\sigma \right]$ denotes the allowable stress. $\left[u\right]=0.234\mathrm{m}$ signifies the admissible value of the greatest displacement of the node located at the top of the jacket. $\left[f_{\max }\right]$ and $\left[f_{\min }\right]$ is the upper limit and lower limit value of the natural frequency, respectively. $u_{max}$ and $H_0$ represent the maximum displacement and natural frequency obtained from the finite element analysis results, respectively. The calculation of stress $\sigma$ is described as follows:

$$\sigma =\sqrt{{\sigma }_a^2+{\sigma }_h^2-{\sigma }_a{\sigma }_h+3{\tau }^2}$$

(52)

where ${\sigma }_h$ denotes the circumferential stress, ${\sigma }_a$ denotes the axial stress, $\tau$ denotes the shear stress.

Furthermore, taking into account the material uncertainties and load uncertainties associated with offshore wind turbine support structures, it is assumed that each variable follows a normal distribution. The uncertain variables and their corresponding properties are illustrated in Table 9.

Table 9.

Uncertain variables.

Random variable	Mean	Cov
Elastic modulus of steel, $E_s\left(\mathrm{GPa}\right)$	210	0.1
Elastic modulus of concrete, $E_c\left(\mathrm{GPa}\right)$	70	0.2
Density of steel, ${\rho }_s(\mathrm{kg}/{\mathrm{m}}^3)$	7850	0.01
Density of concrete, ${\rho }_c(\mathrm{kg}/{\mathrm{m}}^3)$	200	0.05
Buckling strength of steel, $f_y\left(\mathrm{MPa}\right)$	355	0.02
Thrust force, $F_t$ (kN)	781	0.1
Tilting moment, $M_{\mathrm{ti}}$ (kN·m)	38,567	0.1
Torsional moment, $M_{\mathrm{to}}$ (kN·m)	7876	0.1
The mass of tower, $M{s}_t$ (kg)	347,460	0.1

Due to the extremely time-consuming nature of RBDO based on finite element analysis, this numerical example conducts a comparative study on different surrogate model-assisted RBDO methods. Additionally, the RBDO (DO-based) results are also compared. The computational outcomes of the various methods are presented in Table 10.

Table 10.

Comparison of optimization results from different methods.

Method	Number of sample points	Optimal design point	Optimal objective function	LOO error	Iteration count
DO	–	0.966, 0.352, 0.039, 0.045, 0.041, 0.054, 0.062, 0.019, 0.016	65.4719	–	20
PCE	131	1.095, 0.324, 0.039, 0.042, 0.048, 0.054, 0.064, 0.019, 0.020	66.9319	2.5311 × 10−2	41
SPCE	106	1.117, 0.328, 0.041, 0.042, 0.047, 0.055, 0.064, 0.019, 0.020	67.2683	2.0911 × 10−2	37
Kriging	117	1.114, 0.325, 0.040, 0.043, 0.045, 0.053, 0.064, 0.019, 0.020	66.3497	3.1297 × 10−2	39
HDRA	73	1.111, 0.333, 0.039, 0.043, 0.045, 0.054, 0.063, 0.019, 0.020	66.5478	1.0297 × 10−2	26
L-PCK	68	1.127, 0.329, 0.039, 0.043, 0.045, 0.054, 0.063, 0.019, 0.020	66.4374	9.8171 × 10−3	24

As evidenced by the computational results presented in the table, the proposed L-PCK model incurs the least computational cost compared to other surrogate models during its establishment, while simultaneously ensuring the accuracy of the computational results and achieving the smallest error. Compared with HDRA, L-PCK reduced the number of finite-element analyses by ∼7.38%. As illustrated in the table, the error of the L-PCK model is reduced by 68.63% compared to that of the active Kriging model.

To verify the accuracy of the proposed optimization results, finite element analysis was conducted on the design variables obtained from three different optimization designs to assess their constraint violation degrees. As illustrated in Figure 13, the stress distribution contours and deformation contours derived from finite element analysis under various design variables are presented, respectively.

Figure 13.

(a) Presents the finite element analysis verification results of the optimal solution obtained by deterministic optimization; (b) Presents the finite element analysis verification results of the optimal solution obtained by Kriging model-assisted uncertainty optimization; and (c) Presents the finite element analysis verification results of the optimal solution obtained by the L-PCK optimization method proposed in this study.

As illustrated in the figure, the design results obtained through the reliability optimization method exhibit smaller deformations and lower maximum stresses compared to those derived from the deterministic optimization method. Among all the optimized designs, the offshore wind turbine support structure optimized through the L-PCK approach exhibits the minimal deformation and the lowest maximum stress, which substantiates the efficacy and superiority of the L-PCK method.

Conclusion

This study proposes an advanced surrogate modeling strategy assisted by manifold learning and a novel RBDO framework based on this approach. First, the LLE strategy is employed to construct an advanced surrogate model, reducing computational costs and enhancing applicability to complex engineering problems. The PC-Kriging model is utilized to facilitate limit state evaluation during the RBDO process, with the final RBDO framework established based on the PMA. The effectiveness of the proposed RBDO framework is validated through three numerical examples and a case study of an offshore wind turbine support structure, leading to the following conclusions:

1. After implementing active dimensionality reduction via the LLE strategy, the new RBDO framework demonstrates significant time cost reduction compared to conventional RBDO methods when addressing complex engineering problems. In the offshore wind turbine engineering example, the proposed method reduced computational time by 68.83% relative to traditional approaches.

2. The adopted PC-Kriging model exhibits superior computational efficiency and accuracy compared to original single-model surrogate techniques, effectively enhancing the performance of the RBDO framework.

3. Across all mathematical and engineering examples, the proposed RBDO framework consistently demonstrates lower time costs while maintaining accuracy requirements. Relative to RBDO formulations that employ alternative dimension-reduction techniques, the combined LLE–PC-Kriging strategy exhibits markedly higher computational efficiency. This advantage is particularly evident in the final engineering case, confirming the method’s applicability and effectiveness for complex engineering problems.

The research findings indicate that the proposed RBDO method achieves favorable accuracy with superior computational efficiency. However, its generalization capability requires further validation through engineering applications. Future work will focus on exploring more advanced surrogate models, such as Physics-Informed Neural Networks, and investigating optimization design schemes for complex engineering problems involving sophisticated nonlinear constraints or higher-dimensional design spaces.