## Abstract

Deterministic optimization may lead to unreliable design results if significant uncertainty exists. Including reliability constraints in reliability-based design (RBD) can solve such a problem. It is difficult to use current RBD methods to deal with time- and space-dependent reliability when responses vary randomly with respect to time and space. This study employs an envelope method for time- and space-dependent reliability for the optimal design. To achieve high accuracy, we propose an inverse envelope method that converts a time- and space-dependent limit-state function into a time- and space-independent counterpart and then use the second-order saddlepoint approximation to compute the probability of failure. The strategy is to find an equivalent most probable point for a given permitted probability of failure for each reliability constraint. To achieve high efficiency, we use a sequential optimization process to decouple the double-loop structure of RBD. The overall optimization is performed with a sequence of cycles consisting of deterministic optimization and reliability analysis. The constraints of the deterministic optimization are formulated using the equivalent most probable points. The accuracy and efficiency are demonstrated with four examples, including one mathematical problem and three engineering problems.

## 1 Introduction

Engineers always encounter uncertainty in material properties, part dimensions, manufacturing processes, and operational environments [1–3] in all stages of product design and development. Reliability-based design (RBD) is a typical methodology to manage uncertainty by identifying optimal design variables and ensuring satisfied reliability in the design stage [4–6].

RBD minimizes a cost-type objective function while satisfying reliability constraints. If responses are static, meaning that they are time and space invariant, we have static RBD problems, for which many mature methodologies are available, such as double-loop methods [7,8], single-loop single variable approaches [9], sequential optimization method [5,6], and safety factor approaches [10,11]. Du and Chen performed a RBD method by sequential optimization and reliability analysis so that the search of design variables and reliability analysis are executed with a serial of cycles of deterministic optimization and reliability analysis. This helps reduce the computational time [5]. Liang et al. proposed a computationally efficient RBD approach using a single-loop process where the search of the optimum design variables and the reliability analysis are performed simultaneously [9]. Yin and Du developed a modified RBD approach to mechanical component design so that the traditional safety factor design can be used without optimization and complex reliability analysis [10]. Kharmanda demonstrated that the safety factor-based RBD approach is efficient and robust with a new concept of the sequential loop procedure [11]. Tu et al. use a performance measure approach to main robustness and efficiency for evaluating reliability constraints [12].

Many responses are also time-dependent due to time varying stochastic operation conditions and product aging [13]. For instance, the function generator mechanism [14] involves time-dependent motion output. Static RBD methods are not able to handle time-dependent problems. They are extended to time-dependent RBD, and new time-dependent RBD methods have been investigated. In Refs. [15,16], a nested extreme response surface approach accurately carries out time-dependent reliability analysis and determines the optimal designs with good efficiency. The sequential optimization and reliability analysis is extended to time-dependent problems with both stationary stochastic process loads and random loads, and it effectively solves design optimization with dependent reliability constraints [17]. The equivalent most probable point method is proposed to transform the original time-variant RBD problem into an equivalent time-invariant RBD problem formulated by the performance measure approach [18]. The time-dependent concurrent reliability-based design optimization method is developed to improve the confidence of design results with reduced experimental cost and increased efficiency [19]. In Ref. [20], a sequential Kriging modeling approach is introduced to deal with the reliability constraint for time-variant RBD involving stochastic processes.

Although the time-dependent RBD methods have been developed, the most general problems should be addressed where the limit-state functions may take input of stochastic processes, random fields, and tempo-spatial variables. It is still a challenge to address the time- and space-dependent RBD as the research about this topic is limited. Several time- and space-dependent reliability analysis methods have been proposed. Hu and Mahadevan developed a surrogate modeling approach for reliability analysis of a multidisciplinary system [21]. Shi et al. presented a method for the moment estimation of the extreme response using two strategies. One strategy is combining the sparse grid technique and the fourth-moment method while the other is combining the dimensional reduction with the maximum entropy method [22]. Shi and Lu proposed an active learning Kriging method for dealing with dynamic reliability analysis for structure with temporal and spatial parameters [23]. Wei and Du combined first-order reliability method (FORM) and second-order reliability method (SORM) for the time- and space-dependent reliability analysis [24,25]. Yu and Wang developed a general decoupling approach with simulation-based method addressing reliability assessment for time- and space-variant system reliability-based design optimization [26]. Wu and Du extended the envelope method into time- and space-dependent reliability analysis [27,28]. Using the same strategy of the envelope method, we integrate the time- and space-dependent reliability analysis method with optimization in this study [29].

This study extends time- and space-dependent reliability analysis to reliability-based design. The main contributions cover both aspects of analysis and design. (1) Analysis: The new method requires an inverse reliability analysis, whose task is to find an equivalent most probable point (MPP) for a given design and reliability requirement. This is much more challenging than the current direct reliability analysis since the inverse analysis needs to call the direct analysis repeatedly. An efficient analysis method is developed that reduces the number of calls for the direct analysis. (2) Design: We integrate the inverse analysis and optimization by a more efficient means. The new design method performs deterministic optimization and analysis sequentially without embedding the later into the former. During the interactive process, the analysis algorithm performs only partial inverse reliability analysis and converges to the solution of the full inverse analysis at the end of the process. The efficiency is therefore much higher than the direct integration of optimization and inverse analysis.

## 2 Problem Statement

In this section, we define the problem that this study addresses. We also briefly review the sequential RBD, based on which the new method is developed.

### 2.1 Problem Statement.

*i*th response

*y*

_{i}. The associated reliability constraint is that the probability of failure should be smaller than or equal to its allowable value $p~fi$ or $1\u2212R~i$, where [

*R*

_{i}] is the desired reliability.

*n*

_{g}is the number of constraints. In the constraint, the limit-state function is defined by

*m*= 4. (

*z*

_{2},

*z*

_{3}, and

*z*

_{4}represent the

*x*,

*y*, and

*z*coordinates, respectively.) Then, $z=[z1,z2,\u2026,zm]T$ is a vector of the temporal and spatial variables bounded on $\Omega =[z_k,z\xafk]$. Note that the time and spatial variables are quite different physically because they represent different physical aspects. However, they are mathematically the same variables, and it is the reason we present both in the vector $z$. $\Omega $ is a rectangular domain.

*y*is defined by

Solving the time- and space-dependent RBD is time consuming since evaluating *p*_{f} in Eq. (3) is computationally expensive and *p*_{f} must be evaluated repeatedly during the optimization.

### 2.2 Review of Sequential Reliability-Based Design.

Sequential RBD methods in general are more efficient than double-loop methods. It decouples the optimization loop and reliability analysis loop and performs the two loops sequentially. Usually the FORM is used for reliability analysis. FORM searches for the MPP, from which the probability of failure is easily calculated.

The optimal design point is first found from the deterministic optimization loop and then FORM is performed to search for the MPP at this optimal point in the reliability analysis loop. The MPP is then used to modify a reliability constraint for the next deterministic optimization, which is followed by the next reliability analysis. This process repeats until convergence.

*i*th reliability constraint from the reliability analysis. T(·) is the transformation operator for transforming random variables $X$ to independent and standard normal variables $U$. The result of the optimization is the optimal design point $d~$.

*β*

_{i}is the desired reliability index. It is calculated by $\beta i=\u2212\Phi \u22121(p~fi)$, and Φ

^{−1}(·) is inverse cumulative density function.

The final solution can be found after a few cycles of deterministic optimization and reliability analysis. As a result, the efficiency in general is higher than solving the original RBD model directly. Since FORM may not be accurate for highly nonlinear limit-state functions, several studies employ the SORM with higher accuracy and lower efficiency [30,31].

## 3 Sequential Reliability-Based Design With the Envelope Method

### 3.1 Overview.

This gives the solution $u~XMPP=uX$.

The envelope method can be used for reliability analysis when the equivalent MPP $u~XMPP$ is available from Eq. (7). To solve Eq. (7), we at first search for the MPP using FORM at the optimal point of $z$ parameter that minimizes *g*_{i}(·), and then the probability of failure *p*_{fi} is calculated using the envelop method. We update *β*_{i} iteratively until $pfi=p~fi$. The sequential RBD with the envelope method involves cycles of deterministic optimization and equivalent MPP search (reliability analysis).

A design variable could be a distribution parameter of a random variable. For instance, the design task is to determine the length *l* of a shaft, which follows a normal distribution with the mean *μ*_{l} and standard deviation *σ*_{l}. *μ*_{l} is to be determined, and *σ*_{l} is known from the tolerance of *l*. Then we treat *μ*_{l} as a design variable; namely, *d* = *μ*_{l} and *l* = *d* + *X*, where *X* is a random variable that follows a normal distribution with mean of o and standard deviation *σ*_{l}. Design variables are therefore always deterministic in this study.

### 3.2 Time- and Space-Dependent Reliability Analysis.

This subsection discusses the time- and space-dependent reliability using the envelope method. We have already developed a direct envelope method, which calculates the MPP and probability of failure for a given design. The new method in this study is the inverse envelope method, whose task is to find the equivalent MPP for a given probability of failure. The new method is more difficult and time consuming.

There are two cases encountered in the optimization that this study addresses.

*Case 1*: Calling the limit-state function once returns only one response at a specific time instant and spatial point, or a point of $z\u2208\Omega $. Solving for the equivalent MPP in Eq. (5) needs a double-loop procedure. This double procedure will be discussed in Sec. 3.2.1.*Case 2*: Calling the limit-state function once returns all responses at all specific time instants and spatial points, or all points of $z\u2208\Omega $. For instance, if we call a computational fluid dynamics simulation, we obtain all 4D pressure and velocity fields with respect to time and space [3,32]. For this case, solving for the equivalent MPP in Eq. (5) needs only a single loop.

#### 3.2.1 Search for the Equivalent Most Probable Point Using FORM.

There are two constraints in Eq. (7), and directly computing them is too expensive due to a double-loop procedure. We propose a sequential procedure to find the equivalent MPP.

*β*at $z~(0)=[z10,z20,\u2026,zm0]$ with the following model:

Equation (10) produces the initial MPP $uXMPP(1)$. The next analysis is to find the optimal time and space parameter $z~(1)$ by fixing the random variables at its realization $uXMPP(1)$. The next optimal $z~(1)$ can be obtained by another optimization $z~(1)=argminz\u2208\Omega g(d,T(uXMPP(1)),z)$. The details of solving $z~(1)$ are illustrated in Sec. 3.2.3. The process repeats until convergence. This sequential procedure produces the worst-case MPP $uXMPP$ and associated optimal $z~$. After the worst-case MPP is found, the probability of failure with FORM is estimated by $pf=Pr{g(d,T(u),z)\u22640,\u2203z\u2208\Omega}=\Phi (\u2212\Vert uXMPP\Vert )$. However, FORM may not be accurate enough. As a result, we use second-order envelope method to achieve high accuracy.

#### 3.2.2 The Envelope Method.

*g*with respect to

*z*

_{i}.

*U*

_{j}at $uXMPP$ is

*U*

_{j}yields

We use the forward finite difference method with step size $\delta =max(|u|/1000,\epsilon )$, where $\epsilon =10\u22124$, to calculate the derivations in Eq. (22).

The worst-case MPP $uMPP$, gradient $\u2207G$, and the Hessian matrix are now available. The SOSPA is then used to estimate the probability of failure *p*_{fi} after the envelope function is approximated. SOSPA is in general more accurate than FORM because it yields an accurate probability estimation especially in the tail area of a distribution. The details of the implementation of SOSPA are given in Refs. [31,33]. If *p*_{fi} is not equal to $p~fi$, we should update the reliability index *β*_{i} [34], and the MPP search is executed again using Eq. (10). The process is repeated until $pfi=p~fi$. When the probability of failure is equal to the required one, it will produce the equivalent MPP $u~XMPP$. The detailed procedure is given as follows.

*Step 1*: Set the initial reliability index*β*^{(h)}=*β*^{(1)}, the initial optimal point $z(h\u22121)=z(0)$, and the initial MPP $uXMPP(h\u22121)=u0$. The iteration counter*h*= 1.*Step 2*: Search for the inverse MPP using Eq. (10), and obtain the MPP $uMPP(h)$ by solving the following optimization model:(23)${minug(d,T(u),z)s.t.\Vert u\Vert =\beta (h)$*Step 3*: Determine the optimal point $z~(h)$. The optimal point $z~(h)$ minimizes the limit-state function.(24)$z~(h)=argminz\u2208\Omega g(d,T(uXMPP(h)),z)$The optimization method we use in this study is global efficient optimization (EGO) [35].

*Step 4*: Check the convergence criterion.(25)$\Vert uXMPP(h)\u2212uXMPP(h\u22121)\Vert \u2264\epsilon 1$The tolerance $\epsilon 1$ takes a small positive value, for example, 10

^{−4}. If $\Vert uXMPP(h)\u2212uXMPP(h\u22121)\Vert \u2264\epsilon 1$, terminate the iteration. Otherwise, set*h*=*h*+ 1, and return to step 2.*Step 5*: Calculate the gradient $\u2207G$ and Hessian matrix $H$ of the envelope function at $uXMPP(h)$. Calculate the probability of failure*p*_{f}using SOSPA from the above information $uXMPP(h)$, gradient $\u2207G$, and Hessian matrix $H$.

#### 3.2.3 Find the Optimal $z~$.

*z*

_{i}at MPP are obtained as follows:

The optimal point $z~=[z~1,z~2]$ can be obtained by solving Eq. (27). For an explicit and simple limit-state function, the solution of the derivative equations can be obtained analytically. For an implicit and complicated limit-state function, EGO can be used. EGO has been widely used in various areas because it can efficiently search for the global optimum [25].

*h*is the number of initial training points, and the training points of the output dataset are $yin=[g(T(uXMPP),z(i))]i=1,2,\u2026,h$. Once the training dataset $(zin,yin)$ is ready, the next step is to train the initial model using the Gaussian process regression. The initial surrogate model is $y^=g(z)=g(T(uXMPP),z)=F(z)T\gamma +e(z)$, where $F(z)T\gamma $ is a deterministic term, $e(z)$ is a vector of regression functions, $\gamma $ is a vector of regression coefficients, and $e(z)$ is a stationary Gaussian process with zero mean and a covariance is $Cov(e(zi),e(zj))=\sigma e2C(zi,zj)$, where $\sigma e2$ is process variance and

*C*(·, ·) is the correlation function [36]. The initial model may not be accurate; hence new training points are then added one by one so that the model is continuously refined. EGO selects a new training point $znew$ using the expected improvement (EI) metric defined by

*ϕ*(·) is the probability density function.

### 3.3 Extension to Case 2.

*r*is the dimension of $\xi $. Then, the limit-state function becomes $y=g(X~,z)$, where $X~=(X,\xi )$. Thus, the proposed method can still work. Take EOLE as an example for a two-dimensional random field $F(z)$ with $z\u2208(z1,z2)$.

*z*

_{1}and

*z*

_{2}are discretized into $nz1nz2$ points, and the autocorrelation coefficient matrix is given by

*ξ*

_{k}(

*k*= 1, 2, …,

*r*) are independent standard normal variables, $\lambda =(\lambda 1,\lambda 2,\u2026,\lambda r)$ is the eigenvalue vector, and

**ϕ**

_{1},

**ϕ**

_{2}, …,

**ϕ**

_{r}are the corresponding eigenvectors obtained from autocorrelation coefficient matrix $\Sigma $. Note that

*r*is determined as the smallest integer that meets the following criterion:

*η*is a hyperparameter determining the accuracy of the expansion. It takes a value close to, but not larger than 1. The smaller is

*η*, the less accurate is the expansion. If

*η*= 1, the expansion is exact at the points of discretization.

*η*could be set to 0.9–0.99 [37].

### 3.4 Procedure of Sequential Reliability-Based Design With Envelope Method.

After discussing all the details, we now summarize the procedure of the proposed method.

*Step 1*: Initialize parameters.*Step 2*: Perform deterministic optimization. For*k*= 1, solve deterministic optimization at means of random variable. For*k*> 1, perform the following deterministic optimization using the equivalent MPP $u~i,XMPP(k\u22121)$ obtained from the (*k*− 1)th cycle. The solution is $d(k)$.(35)${mind,\mu Xf(d)s.t.gi(d,T(ui,XMPP(k\u22121)),z)\u22640,i=1,2,\u2026,ndL\u2264d\u2264dU$*Step 3*: Perform time- and space-dependent reliability analysis at $d(k)$ for each constraint. At first, search for the equivalent MPP given*β*^{(k)}. Obtain $u~XMPP(k)$, gradient $\u2207G(u~XMPP(k))$, and Hessian matrix $Hu~XMPP(k)$. Note that if the inputs of limit-state function are random variables and random fields, the method in Sec. 3.3 is used to find the $u~XMPP(k)$. Next, calculate the probability of failure $pfi(k)$ using SOSPA.*Step 4*: Check the convergent criterion by(36)$\epsilon =|pfi(k)\u2212[Pfi][Pfi]|\u2264\epsilon 2$If convergence is reached, the optimal solution is found at $d(k)$, and the process stops. Otherwise update the reliability index*β*^{(k+1)}by(37)$u~XMPP(k+1)=\beta (k+1)\u2207g(u~XMPP(k))\Vert u~XMPP(k)\Vert $The flowcharts of the proposed method are given in Figs. 1 and 2.

The new method has the same limitations as other MPP-based reliability methods. Its accuracy will deteriorate if there are multiple MPPs, which occur when the envelope function is piecewise and indifferentiable. We could use the method of composite limit-state function for time-dependent reliability to solve this problem [38]. Its extension to time- and space-dependent reliability needs further investigation.

The overall process can typically converge within five cycles. Since the process starts from deterministic optimization in the first cycle with mean values of random variables, the reliability of an active constraint is approximately 50% at the optimal point, much lower than the required reliability. In the following reliability analysis, the MPP is found for the given design point and the required reliability. Then in the deterministic optimization of the second cycle, the MPP is plugged into the corresponding constraint function, which guarantees that the reliability will be improved and much closer to the target reliability, but still smaller than the target reliability. Then the reliability analysis is conducted again. More cycles are performed until convergence. Hence the design process always starts from the infeasible region, and the design point will gradually approach constraint boundaries and the final optimal point. On other words, the reliability will gradually be improved until reaching the target. The entire design process will likely converge if both deterministic optimization and reliability analysis converge. The method, however, may not work if the reliability analysis does not identify a global MPP; for instance, a local MPP is found in one and a global MPP is found in another cycle.

## 4 Examples

*p*

_{f}is the result from a non-MCS method while $pfMCS$ is from MCS.

Since case 1 is more complicated by case 2 and case 1 is the focus of this study, all the examples discussed here belong to case 1.

### 4.1 Example 1: A Mathematical Problem.

*X*

_{1}and

*X*

_{2}are normally distributed with $X1\u223cN(\mu X1,0.6)$ and $X2\u223cN(\mu X2,0.6)$. The time

*t*varies over the interval [0,1], and the spatial variable

*s*changes over the interval [7]. The design variables are $\mu X1$ and $\mu X2$. The limit-state function is defined by

The allowable reliability index *β* is 3. The optimal results are achieved by following the detailed steps of the proposed method. The proposed method produces the following results: design variables are *w* = 2.7177 and *h* = 4.5691; the objective function is *f* = −7.2867. The results by SORM/DL and FORM/DL are provided in Table 1. The proposed method and SORM/DL have the same accuracy as their errors are the same. This shows that the former method converges to the same solution from the latter method; the former method is much more efficient than the latter method as indicated by the number of function calls *N*_{call} in Table 1. As expected, the proposed method is more accurate than FORM/DL that uses a first-order approximation. The proposed method is also more efficient than FORM/DL although the former and latter methods use second and first approximations, respectively. The convergence history is shown in Fig. 3.

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. | −7.2867 | −7.2867 | −7.3096 |

μ | (2.7177, 4.5691) | (2.7177, 4.5691) | (2.7275, 4.5821) |

$pfMCS$ (10^{−3}) | 1.3679 | 1.3666 | 1.4942 |

Error (%) | 1.33 | 1.23 | 10.68 |

N_{calls} | 341 | 4317 | 747 |

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. | −7.2867 | −7.2867 | −7.3096 |

μ | (2.7177, 4.5691) | (2.7177, 4.5691) | (2.7275, 4.5821) |

$pfMCS$ (10^{−3}) | 1.3679 | 1.3666 | 1.4942 |

Error (%) | 1.33 | 1.23 | 10.68 |

N_{calls} | 341 | 4317 | 747 |

### 4.2 Example 2: Three Bar System.

A truss structure is shown in Fig. 4. Each bar of the system has its cross-sectional area *A*_{i} and modulus of elasticity, *E*_{i}, *i* = 1, 2, 3. The coefficient of thermal expansion of all bars is $\alpha =12\xd710\u22126\u2218C\u22121$. The temperature change is related to the installation height of the truss structure and is given by $\Delta T=Te\u22120.01(\Delta h2+2\Delta h\u22122)2$, where Δ*h* ∈ [1, 6] m is the difference between two different installation heights. A downward force *P* = *P*_{0}(0.9 + 0.1cos(0.2*t*)) is applied at joint *A*, where *t* ∈ [0, 10] years. All the random variables are given in Table 2. The design variables are the cross-sectional areas of the bars $\mu A1$, $\mu A2$, and $\mu A3$.

Variable | Mean | Standard deviation | Distribution |
---|---|---|---|

A_{1}(mm^{2}) | $\mu A1$ | 0.6 | Normal |

A_{2}(mm^{2}) | 0.6 | Normal | |

A_{3}(mm^{2}) | $\mu A3$ | 0.6 | Normal |

E_{1}(GPa) | 200 | 20 | Normal |

E_{2}(GPa) | 200 | 20 | Normal |

E_{3}(GPa) | 200 | 20 | Normal |

P_{0}(KN) | 40 | 6 | Normal |

L_{AD}(mm) | 200 | 2 | Normal |

L_{AB}(mm) | 231 | 2.31 | Normal |

L_{AC}(mm) | 283 | 2.83 | Normal |

$T(\u2218C)$ | 35 | 7 | Normal |

σ_{yield}(GPa) | 7.5 × 10^{8} | 4 × 10^{7} | Normal |

Variable | Mean | Standard deviation | Distribution |
---|---|---|---|

A_{1}(mm^{2}) | $\mu A1$ | 0.6 | Normal |

A_{2}(mm^{2}) | 0.6 | Normal | |

A_{3}(mm^{2}) | $\mu A3$ | 0.6 | Normal |

E_{1}(GPa) | 200 | 20 | Normal |

E_{2}(GPa) | 200 | 20 | Normal |

E_{3}(GPa) | 200 | 20 | Normal |

P_{0}(KN) | 40 | 6 | Normal |

L_{AD}(mm) | 200 | 2 | Normal |

L_{AB}(mm) | 231 | 2.31 | Normal |

L_{AC}(mm) | 283 | 2.83 | Normal |

$T(\u2218C)$ | 35 | 7 | Normal |

σ_{yield}(GPa) | 7.5 × 10^{8} | 4 × 10^{7} | Normal |

*A*, denoted as Δ

*δ*, is greater than the allowable displacement

*δ*

_{0}, and the limit-state function is defined by

*δ*

_{0}= 0.64, and the perpendicular displacement of joint

*A*is calculated by

The allowable reliability indexes are *β*_{i} = 2.3, *i* = 1, 2. The results from three are given in Table 3. The solution of the proposed method is ** μ** = [53.41, 38.83, 60] mm

^{2}, which satisfy all the design constraints. The probabilities of failure at the optimal design points are $pf1MCS=0.0107$ and $pf2MCS=0.0107$, which is the same as the probabilities of failure with SORM/DL method. Besides, the proposed method has higher efficiency than SORM/DL in terms of the

*N*

_{calls}. Consequently, the proposed method achieves the most accurate and efficient results.

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. ( × 10^{4} mm^{3}) | 3.6616 | 3.6624 | 3.7189 |

μ (mm) | (53.41, 38.83,60) | (53.31, 38.95,60) | (53.05, 38.57,60) |

$pf1MCS$ | 0.0107 | 0.0107 | 0.0148 |

$\epsilon pf1$ (%) | 0.25 | 0.46 | 4.89 |

$pf2MCS$ | 0.0107 | 0.0107 | 0.0112 |

$\epsilon pf2$ (%) | 0.24 | 0.46 | 4.89 |

N_{calls} | 2194 | 23,500 | 7264 |

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. ( × 10^{4} mm^{3}) | 3.6616 | 3.6624 | 3.7189 |

μ (mm) | (53.41, 38.83,60) | (53.31, 38.95,60) | (53.05, 38.57,60) |

$pf1MCS$ | 0.0107 | 0.0107 | 0.0148 |

$\epsilon pf1$ (%) | 0.25 | 0.46 | 4.89 |

$pf2MCS$ | 0.0107 | 0.0107 | 0.0112 |

$\epsilon pf2$ (%) | 0.24 | 0.46 | 4.89 |

N_{calls} | 2194 | 23,500 | 7264 |

### 4.3 Example 3: A Cantilever Beam.

*F*

_{1}and

*F*

_{2}. The length of the cantilever beam

*L*is 100 in. The objective is to minimize the volume

*f*=

*μ*

_{w}

*μ*

_{h}

*L*, where

*w*and

*h*represent the width and height of the beam cross section, respectively. There are two failure modes. The first failure mode is that the stress at the fixed end is greater than the allowable yield stress

*S*

_{y}, and the second failure mode is that the tip displacement of the beam is greater than the allowable displacement

*D*

_{0}= 2.5 in. The two limit-state functions are given by

*F*

_{1}is a non-stationary Gaussian random field, whose mean is $\mu F1=500e0.01((s\u2212(1/2))2+(t\u22126)2)lbs$ and standard deviation is $\sigma F1=50lbs$. The spatial variable is

*s*∈ [0, 1] in. and temporal variable is

*t*∈ [0, 10] years. The autocorrelation coefficient function is given by

Variable (unit) | Mean | Standard deviation | Distribution |
---|---|---|---|

w(in.) | μ_{w} | 0.01 | Normal |

h(in.) | μ_{h} | 0.01 | Normal |

F_{2}(lb) | 1 × 10^{3} | 1 × 10^{2} | Normal |

E(psi) | 2.9 × 10^{7} | 1 × 10^{5} | Normal |

S_{y}(psi) | 3.9 × 10^{4} | 500 | Weibull |

Variable (unit) | Mean | Standard deviation | Distribution |
---|---|---|---|

w(in.) | μ_{w} | 0.01 | Normal |

h(in.) | μ_{h} | 0.01 | Normal |

F_{2}(lb) | 1 × 10^{3} | 1 × 10^{2} | Normal |

E(psi) | 2.9 × 10^{7} | 1 × 10^{5} | Normal |

S_{y}(psi) | 3.9 × 10^{4} | 500 | Weibull |

The allowable reliability indexes are *β*_{i} = 3, *i* = 1, 2. Table 5 shows that the optimal design variables are *w* = 3.9541 in. and *h* = 2.2531 in. in, and the objective function is *f* = 890.9152 in.^{3} by the proposed method. The probabilities of failure obtained at the optimal design variables by MCS are $pf1MCS=1.3633\xd710\u22124$ and $pf2MCS=1.3645\xd710\u22124$. The results are more accurate than those of SORM/DL and FORM/DL methods. The proposed method is the most efficient method as the number of function calls is 1358 compared with those of SORM/DL and FORM/DL, which are 12,904 and 2098, respectively. In general, the proposed method is the best method for this example due to its high accuracy and efficiency.

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. (in.^{3}) | 889.3705 | 889.3766 | 890.0996 |

μ (in.) | (3.9375, 2.2587) | (3.9373, 2.2589) | (3.9388, 2.2573) |

$pf1MCS(\xd710\u22124)$ | 1.3435 | 1.3429 | 1.7410 |

$\epsilon pf1$ (%) | 0.4740 | 0.5184 | 0.918 |

$pf2MCS$ ( × 10^{−4}) | 1.3399 | 1.3348 | 1.4195 |

$\epsilon pf2$ (%) | 0.7407 | 1.1185 | 5.1561 |

N_{calls} | 1358 | 12,904 | 2098 |

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. (in.^{3}) | 889.3705 | 889.3766 | 890.0996 |

μ (in.) | (3.9375, 2.2587) | (3.9373, 2.2589) | (3.9388, 2.2573) |

$pf1MCS(\xd710\u22124)$ | 1.3435 | 1.3429 | 1.7410 |

$\epsilon pf1$ (%) | 0.4740 | 0.5184 | 0.918 |

$pf2MCS$ ( × 10^{−4}) | 1.3399 | 1.3348 | 1.4195 |

$\epsilon pf2$ (%) | 0.7407 | 1.1185 | 5.1561 |

N_{calls} | 1358 | 12,904 | 2098 |

### 4.4 Example 4: Thermal Deflection of a Bimetallic Beam.

Thermal expansion or contraction of a bimetallic beam occurs due to temperature change. The temperature in the operating room varies during the day and night and is given by Δ*T* = *T*(0.01sin(0.1*t*) + 1), where *t* ∈ [0, 24] h. The typical bimetallic beam consists of two materials bonded together: copper and invar. *E*_{C} is Young’s modulus of copper, and *E*_{I} is Young’s modulus of invar. The length of the beam depends on the location of the installation, which is given by *L* = *L*_{0}(− *s*^{2} + *s* + 1), where *s* ∈ [0, 1] m. When the temperature change as a thermal load applies on the beam, the beam will deflect in the perpendicular direction at the right end side shown in Fig. 6. The design variables are $d=(h,w)$, where *h* and *w* represent the height and width of the cross-sectional area of the beam, and their means are *μ*_{h} and *μ*_{w}, which are to be determined. All the random variables are listed in Table 6.

Variable | Mean | Standard deviation | Distribution |
---|---|---|---|

w(m) | μ_{w} | 5 × 10^{−4} | Normal |

h(m) | μ_{h} | 5 × 10^{−4} | Normal |

L(m) | 1 × 10^{−1} | 1 × 10^{−3} | Normal |

E_{C}(Pa) | 1.37 × 10^{11} | 1.37 × 10^{7} | Lognormal |

E_{I}(Pa) | 1.30 × 10^{11} | 1.3 × 10^{7} | Lognormal |

$T(\u2218C)$ | 130 | 13 | Lognormal |

Variable | Mean | Standard deviation | Distribution |
---|---|---|---|

w(m) | μ_{w} | 5 × 10^{−4} | Normal |

h(m) | μ_{h} | 5 × 10^{−4} | Normal |

L(m) | 1 × 10^{−1} | 1 × 10^{−3} | Normal |

E_{C}(Pa) | 1.37 × 10^{11} | 1.37 × 10^{7} | Lognormal |

E_{I}(Pa) | 1.30 × 10^{11} | 1.3 × 10^{7} | Lognormal |

$T(\u2218C)$ | 130 | 13 | Lognormal |

*δ*= 8 × 10

^{−3}. The limit-state function is given by

**d**,

*E*

_{C},

*E*

_{I}, Δ

*T*) is solved by the finite element method, and matlab partial differential equation toolbox is used to solve this example.

All the results are listed in Table 7. The proposed method is the most accurate method among the three methods as the error is only 3.04% as it uses second-order saddlepoint to find the optimal design results. It is the most efficient method as well with the evidence that the number of function calls is only 332 since it uses the sequential loop to improve the computation efficiency.

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. ( × 10^{−6}m^{2}) | 1.9197 | 1.9181 | 1.9312 |

μ (10^{−4}) | (8, 2.3996) | (8, 2.3975) | (8, 2.4319) |

$pfMCS$ | 0.0013 | 0.0014 | 0.0012 |

Error (%) | 3.04 | 5.71 | 10.43 |

N_{calls} | 332 | 3645 | 573 |

Method | New method | SORM/DL | FORM/DL |
---|---|---|---|

Obj. ( × 10^{−6}m^{2}) | 1.9197 | 1.9181 | 1.9312 |

μ (10^{−4}) | (8, 2.3996) | (8, 2.3975) | (8, 2.4319) |

$pfMCS$ | 0.0013 | 0.0014 | 0.0012 |

Error (%) | 3.04 | 5.71 | 10.43 |

N_{calls} | 332 | 3645 | 573 |

## 5 Conclusions

This study develops a new sequential RBD with the envelope method for time- and space-dependent reliability. The challenge in this work is to search for the equivalent most probable point (MPP), which can be found by iterating the MPP search and updating the equivalent reliability index. When limit-state functions have input with random variables and the temporal/spatial domain Ω, a single-loop MPP search is used which is much more efficient than the sequential loop approach. Once the equivalent MPP is available, the time- and space-dependent problem is transformed into a static counterpart and the second-order saddlepoint approximation is used to estimate the reliability with higher accuracy. The equivalent MPP assures that the overall optimization is performed sequentially in cycles of deterministic optimization and reliability analysis. The proposed strategy has been proven to be effective in four examples.

The proposed method is more accurate than the first-order methods since it uses the second-order saddlepoint approximation to estimate the reliability. The new method, however, shares the same limitations as the MPP-based methods. Its accuracy will deteriorate if multiple MPPs exist. The other limitation is that the method cannot handle the case where the MPP occurs on boundaries of the time and space domain. In this case, the accuracy of reliability prediction will also deteriorate. The future work is to use a system reliability method to accommodate multiple MPPs and derive a new algorithm for an MPP on the boundaries of the time and space domain.

## Acknowledgment

The support from National Science Foundation under grant under Grant No. 1923799 is acknowledged.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.