2.1. System reliability with independent component states

A series system is shown in Figure 1, in which the components in the system are denoted by C ₁, C ₂, … , C _n . The system will fail if one of its components fails. If all the component failures are independent, the system reliability R _S is

(1) $$R_{\rm S} = \prod\limits_{i = 1}^n {R_i}\comma $$

where R _i (i = 1, 2, … , n) is the reliability of component i.

Fig. 1. A series system.

Component reliability can be estimated by a statistics-based approach with testing or field data. It can also be estimated by a physics-based approach. If the latter approach is used, component reliability is given by

(2) $$R = \Pr \left\{ {Y = g({\bf X}) \gt 0} \right\}\comma \,$$

where X is a vector of random input variables, and Y is the state variable. If Y > 0, the component functions; otherwise, the component fails.

In this work, we focus on mechanical applications where series systems are usually involved.

2.2. System reliability bounds

Eq. (1) is easy to use, but may produce a large error due to the independent component assumption and may be too conservative. The actual system reliability is bounded as shown (Arnljot & Rausand, Reference Arnljot and Rausand2009):

(3) $$\prod\limits_{i = 1}^n {R_i \le R_{\rm S} \le \min \{ R_i \}}\comma \quad i = 1\comma \, \ldots\comma \, n. \,$$

If a mechanical system consists of 20 components with identical component reliabilities R = 0.999, Eq. (3) gives the bounds of 0.9802 ≤ R _S ≤ 0.999. The bounds may be too wide to help system designers to compare design concepts for concept selection.

2.3. System reliability with components sharing the same system load

To improve system reliability analysis, we performed a preliminary study for systems whose components share the same stochastic system load L (Hu & Du, Reference Hu and Du2016). The system designers have good knowledge about L and therefore know the cumulative distribution function (CDF) of L; L is distributed through components, and the component load L _i (i = 1, 2, … , n) of component i is a function of L. Such a function is assumed to be

(4) $$L_i = w_i L\comma \,$$

where w _i indicates the fraction of the load that the component shares. This w _i can be determined from a system level analysis, such as a force analysis.

System designers request component designers to provide component reliability functions at different component load levels, specified by variable l. The component designers may conduct experiments or use a physics-based approach to calculate component reliability R _i by varying the values of l. Then the component reliability functions R _i (l) are available to system designers.

System designers then assume that the component state could be predicted by the following component limit-state function:

(5) $$g_i (L_i ) = Y_i = S_i - L_i = S_i - w_i L\comma \,$$

where S _i is the general resistance of the component. The component limit-state function should reproduce the same component reliability, namely,

(6) $$R_i (L_i ) = \Pr \{ g_i (L_i ) \gt 0\}. $$

The probability of system failure is then given by

(7) $$\eqalign{ p_{\rm fs} &= {\rm Pr} \{S_1 \lt L_1 \bigcup S_2 \lt L_2 \bigcup \cdots \bigcup S_n \lt L_n\} \cr &= {\rm Pr} \left\{ \mathop \bigcup \limits_{i = 1}^n S_i \lt w_i L \right\}. } $$

The system reliability is then available and is given by

(8) $$R_{\rm S} = 1 - p_{\,\rm fs}. $$

It is obvious that component failure events S _i < w _i L are dependent because of the common random variable L. The component dependence is therefore considered automatically. The reliability function R _i (l) is directly related to the CDF of S _i , because the CDF of S _i is 1 − R _i (l) (Hu & Du, Reference Hu and Du2016). If S _i and L are independent, system designers know the joint distribution of all the random variables in Eq. (7), and thus, they can use Eq. (7) to find the system reliability.

3.1. Procedure of system reliability prediction

To make system reliability prediction possible, system designers ask component suppliers to provide component reliabilities with respect to their component loads. As the information of component reliabilities may be in different forms, for system designers, the first step is to formulate component reliability functions R _i (l), i = 1, 2, … , n, with respect to the component load l. The second step is to construct composite component limit-state functions g _i (·) based on R _i (l). A composite component limit-state function ensures that it can reproduce accurate component reliability regardless of the number of failure modes the component may have. A composite component limit-state function does not require any component design details, and thus, it prevents the proprietary information of the component supplier. Because the common system load appears in all composite component limit-state functions, the dependence between component states is preserved. This helps improve the accuracy of system reliability prediction. The last step is to perform system reliability analysis. The flowchart indicating the procedure is given in Figure 2.

Fig. 2. Flowchart of the proposed method.

3.2. Formulation of component reliability functions

Component designers may use different methods to estimate component reliabilities R _i (l), such as using testing, field data, simulations, or a physics-based method. This may result in different forms of information about R _i (l), such as limited reliability data, a scatter plot, or a mathematical model. If no mathematical model exits, system designers need to fit a model from these limited data. As will be discussed in Section 3.3, the probability of failure, p _fi (l) = 1 − R _i (l), is actually the CDF of the general component resistance. Then, the task becomes to fit a CDF model. Many methods could be used for the CDF fitting such as metamodeling methods, the saddle point approximation, and the Weibull analysis.

Now we discuss how system designers could fit a CDF model given the limited reliability data. For a general component with probability of failure p _f (l), the available data are given as a set of (l _j , p _f(l _j )), j = 1, 2, … , m. Let the continuous mathematical model be

(9) $$p_{\rm f} = H(l).$$

Next, we discuss two specific approaches to obtain H(l) from (l _j , p _f(l _j )), j = 1, 2, … , m.

3.2.1. Kriging method

The kriging method has been widely used in engineering applications, including reliability analysis (Kleijnen, Reference Kleijnen2009; Viana & Haftka, Reference Viana and Haftka2012; Kolios & Salonitis, Reference Kolios and Salonitis2013). The kriging method considers the mathematical model in Eq. (9) as a realization of a Gaussian process given by (Sacks et al., Reference Sacks, Welch, Mitchell and Wynn1989):

(10) $$p_{\rm f} (l) = H(l) = {\rm \alpha} (l){\rm \xi} + Z(l)\comma \,$$

where α(l) is a regression function, ξ is the regression coefficient, and Z(·) is a stationary Gaussian process with zero mean. The covariance between two points l _i and l _j is defined by

(11) $$\eqalign{ {\rm cov}\lsqb Z\lpar {l_i} \rpar \comma \,Z\lpar {l_j} \rpar \rsqb = {\rm \sigma} _Z^2 K(l_i\comma \, &l_j )\comma \, \cr &i = 1\comma \,2\comma \, \ldots\comma \,m; \; j = 1\comma \,2\comma \, \ldots\comma \, m\comma \,}$$

where ${\rm \sigma} _Z^2 $ is the variance of the Gaussian process, and K( · , · ) is the correlation function and is commonly defined by the following Gaussian correlation (Sacks et al., Reference Sacks, Welch, Mitchell and Wynn1989; Lophaven et al., Reference Lophaven, Nielsen and Sondergaard2002):

(12) $$K\lpar {l_i\comma \, l_j} \rpar = \exp \lsqb { - {\rm \theta} \lpar {l_i - l_j} \rpar^2} \rsqb \comma \,$$

where θ is a parameter that indicates the correlation between the points. The best linear unbiased predictor (Sacks et al., Reference Sacks, Welch, Mitchell and Wynn1989) of H(l) gives to a random prediction

(13) $$\hat p_{\rm f} = \hat H(l) \sim N\lpar {{\rm \mu}_H \lpar l \rpar \comma \,{\rm \sigma}_H^2 \lpar l \rpar } \rpar \comma \,$$

where the prediction μ _H (·) and the associated variance are computed by

(14) $${\rm \mu}_H (l) = {\rm \alpha} (l)\hat \xi + {\bf r}^T (l){\bf K}^-1 ({\bf p}_{\bf f} - {\bf F}{\rm \hat \xi} )\comma \,$$

(15) $${\rm \sigma} _H^2 (l) = {\rm \hat \sigma} _Z^2 \left\{ \matrix{1 - [{\bf r}(l)]^T {\bf K}^{ - 1} {\bf r}(l) \hfill \cr + \left[ {{\bf F}^T {\bf K}^{{\rm-1}} {\bf r}{\rm (}l{\rm )} - {\rm \alpha (}l{\rm )}} \right]^T \left( {{\bf F}^T {\bf K}^{{\rm-1}} {\bf F}} \right)^{{\rm-1}} \hfill \cr \left[ {{\bf F}^T {\bf K}^{{\rm-1}} {\bf r}{\rm (}l{\rm )} - {\rm \alpha (}l{\rm )}} \right] \hfill} \right\}\comma \,$$

in which K is the correlation matrix defined by ${\bf K} = \lsqb {K(l_i\comma \, l_j )} \rsqb \comma \, {\bf p}_{\bf f} $ is a column vector of responses of current sample points, and ${\bf r}( \cdot )$ is the vector of cross-correlations between the m samples and the prediction point, ${\bf r}(l) = \lsqb {R(l\comma \,l_1 )\comma \, \ldots\comma \, R(l\comma \,l_m )} \rsqb ^T $ . F is a column vector with rows α(l _i ), i = 1, 2, … , m, and ${\rm \hat \sigma} _Z^2 $ is the maximum likelihood estimation of the process variance,

(16) $${\rm \hat \sigma} _Z^2 = \displaystyle{1 \over m}({\bf p}_{\bf f} - {\bf F}{\rm \hat \xi} )^T {\bf K}^{ - 1} ({\bf p}_{\bf f} - {\bf F}{\rm \hat \xi} )\comma \,$$

and ${\rm \hat \xi} $ is the generalized least square estimate of ξ,

(17) $${\rm \hat \xi} = \left[ {{\bf F}^T {\bf K}^{ - 1} {\bf F}} \right]^{ - 1} {\bf F}^T {\bf K}^{ - 1} {\bf p}_{\bf f}. $$

Substituting Eqs. (14) and (15) into Eq. (13), system designers obtain the reliability mathematical model in the form of $p_{\rm f} = \hat H(l)$ for p _f = H(l) in Eq. (9).

3.2.2. Weibull method

A Weibull distribution can fit different data and distributions. Due to this advantage, system designers may use a Weibull model to fit the component reliability data. A three-parameter Weibull distribution is given by

(18) $$p_{\rm f} = H(l) = 1 - \exp \bigg( { - \bigg( {\displaystyle{{l - {\rm \gamma}} \over {\rm \eta}}} \bigg)^{\rm \beta}} \bigg)\comma \,$$

in which l > γ, β > 0, η > 0. The location parameter γ defines the location of the distribution; β is the shape parameter, and η is the scale parameter.

For a given set of (l _j , p _f(l _j )), j = 1, 2, … , m, system designers could use the maximum likelihood method (Lockhart & Stephens, Reference Lockhart and Stephens1994) to find the three distribution parameters. They could also use a curve fitting method (Hamming, Reference Hamming2012) to find the three distribution parameters.

Other regression analysis methods could also be used for the CDF fitting. One example showing the CDF fitting follows. Suppose the probabilities of component failure p _f at seven load levels are given and are shown in Figure 3. A mathematical model of p _f with respect to the component load l can be then fitted as shown in Figure 4.

Fig. 3. Component reliability data.

Fig. 4. The complete p _f model.

3.3. Reconstruction of component limit-state functions

The next task of system designers is to reconstruct component limit-state functions, which should meet the following requirements:

• • Do not require component design details

• • Maintain dependence between component states

• • Be functions of the system load

• • Be easy to evaluate

• • Accommodate multiple component failure modes

Based on these requirements, for the ith component, system designers reconstruct the limit-state function in the form of

(19) $$Y_i = S_i - L\comma \,$$

where S _i is the general component resistance, and L is the system load. Note that no matter how many failure modes the component may have, there is only one reconstructed component limit-state function as shown in Eq. (19).

Although the reconstructed component limit-state function is linear with respect to L, it can accommodate the situation where the actual component limit-state function is nonlinear with respect to L. One example follows. Let the yield strength of the i-th component be S _y . If the maximum stress is h(L), where h(·) is a nonlinear function, also depending on other component parameters, such as dimensions, and then component designers build their limit-state function as

(20) $$Y^{\prime}_i = S_y - h(L).$$

Theoretically, they can solve for L by letting S _y − h(L) = 0 at the limit state and obtain

(21) $$L = h^{ - 1} (S_y )\comma \,$$

where h ⁻¹(·) is the inverse function of h(·).

Then the limit-state function is modified as

(22) $$Y_i = h^{ - 1} (S_y ) - L.$$

Let S _i = h ⁻¹(S _y ), which is regarded as the general component resistance. Then Eq. (22) is exactly the one reconstructed by system designers in Eq. (19). This indicates that the reconstruct component limit-state functions do cover actual component limit-state functions that are nonlinear with respect to the system load.

Before explaining the procedure of reconstructing the composite limit-state function, we first prove that the probability of component failure p _fi (l) is the CDF of the general component resistance S _i . According to Eq. (19), for a constant l,

(23) $$p_{\,{\rm f}i} (l) = \Pr (S_i \lt l).$$

The CDF of S _i is defined by

(24) $$F_{S_i} (s) = \Pr (S_i \lt s).$$

Replacing l with s in Eq. (23), we have $p_{\,fi} (s) = \Pr (S_i \lt s).$ As a result,

(25) $$F_{S_i} (s) = p_{\,{\rm f}i} (s) = 1 - R_i (s).$$

Because system designers know the component reliability function R _i (l) or probability of component failure p _fi (l), they also know the CDF of the general component resistance S _i .

The composite component limit-state function is not only a simple (linear) function but also safeguards the proprietary information of component designers. Next, let us look at the component design of the example that will be presented in Section 4.

In the example, Component 2 has two failure modes due to excessive normal stress and excessive shear stress. The component designer decides to use a physics-based approach to evaluate the component reliability. The limit-state functions of the two failure modes are given by

(26) $$Y_{21} = S_y - \displaystyle{{(h + H_1 )} \over {W_x}} \displaystyle{L \over 2}\comma \,$$

(27) $$Y_{22} = {\rm \tau} - \displaystyle{1 \over {hb}}\displaystyle{L \over 2}\comma \,$$

in which L/2 is the load shared by the component, and h, H ₁, W _x , and b are random parameters related to component details. The two limit-state functions indicated that component details are required for the component reliability analysis. The details include material properties, component structure, and component dimensions.

The two limit-state functions for the two failure modes can be rewritten as

(28) $$Y^{\prime}_{21} = \displaystyle{{2W_x S_y} \over {h + H_1}} - L = S_{21} - L\comma \,$$

(29) $$Y^{\prime}_{22} = 2{\rm \tau}hb - L = S_{22} - L\comma \,$$

in which $S_{21} = (2W_x S_y)/(h + H_1)$ and S ₂₂ = 2τhb. Then, the probability of component failure is

(30) $$\eqalign{ p_{\,{\rm f}2} &= \Pr \!\lcub {S_{21} \lt L \bigcup S_{22} \lt L} \rcub \cr &= \Pr\! \lcub {\min \lpar {S_{21}\comma \, S_{22}} \rpar \lt L} \rcub = \Pr\! \lcub {S_2 \lt L} \rcub \comma \,}$$

where S ₂ = min(S ₂₁, S ₂₂) is the general component resistance in Eq. (19). Note that the details such as h, H ₁, W _x , and b in Eqs. (26) through (29) are only known to component designers who could find the component probabilities of failure at different load levels of load l _i (i = 1, 2, … , n) by testing or using a physics-based reliability approach. Then the results could be provided to the system designers in the form of (p _f2(l _i ), l _i ). As shown in Eqs. (24) and (25), the probability of component failure p _f2 is exactly the CDF of the general component resistance S ₂. Thus, the distribution of S ₂ is known to system designers, and with this distribution, they no longer need any design details. No proprietary information is therefore required. Equation (19) also indicates that the system load L appears in all the reconstructed component limit-state functions, and the dependence between component states is automatically maintained. Meanwhile, the composite limit-state function takes into account the multiple component failure modes, and it has a simple expression to evaluate. Thus, the obtained composite limit-state functions satisfy all the requirements mentioned above.

3.4. System reliability analysis

We now discuss how system designers use the reconstructed composite component limit-state functions in Eq. (19) to predict system reliability. The probability of system failure is given by

(31) $$p_{\rm fs} = \Pr \left\{ {\bigcup\limits_{i = 1}^n {Y_i = S_i - L \lt 0}} \right\}.$$

The prerequisite for calculating p _fs is to find the joint probability distribution of S _i (i = 1, 2, … , n) and L. We now discuss how to obtain such a joint probability distribution.

Denote the n + 1 input random variables by Z = (S ₁, S ₂, … , S_n , L) and all the output variables by Y = (Y ₁, Y ₂, … , Y_n ). As discussed in Section 3.3, the general component resistances S _i (i = 1, 2, … , n) are determined by component material properties, concrete component structures, geometric dimensions, and other component parameters. As all the components are independently designed, manufactured, and tested by different suppliers, their general resistances are likely statistically independent. The system load L is also independent from the general component resistances. Thus, all the components in Z are independent.

Denote the CDF of S _i and L by $F_{S_i} (s_i )$ , and F _L (l), respectively. The joint CDF of Z is then given by

(32) $$F_{\bf Z} ({\bf z}) = F_L (l)\prod\limits_{i = 1}^n {F_{S_i} (s_i )}\comma \, $$

where ${\bf z} = (s_1\comma \, s_2\comma \, \ldots\comma \, s_n\comma \, l)$ . Because system designers know CDFs of S _i and L, it is easy for them to predict the probability of system failure. Denote the joint probability density function (PDF) of Z by $f_{\bf Z} ({\bf z})$ . Then the probability of system failure is computed by

(33) $$p_{\rm fs} = \int\limits_\Omega {\,f_{\bf Z} ({\bf z})} d{\bf z}\comma $$

where Ω is the system failure region defined by

(34) $$\Omega = \left\{ {{\bf Z} \vert S_i \lt L\comma \quad i = 1\comma \,2\comma \, \ldots\comma \,n} \right\}\comma \,$$

where p _fs can be calculated by an numerical integration or Monte Carlo simulation (MCS). Next, we demonstrate this with two special cases.

In the first case, the components S _i (i = 1, 2, … , n) of Z = (S ₁, S ₂, … , S_n , L) follow Weibull distributions, while L follows a distribution with PDF F _L (l). Note that Y _i = S _i − L, and thus Y ₁, Y ₂, … , Y _n are dependent. Because the distribution parameters of Y are unknown, it is not possible to directly find p _fs using Eq. (31). However, as we know that S ₁, S ₂, … , S _n , and L are independent, the joint PDF of Z could be calculated by

(35) $$\eqalign{ f_{\bf Z} ({\bf z}) &= f_{\bf Z} (s_1\comma \,s_2\comma \, \ldots\comma \,s_n\comma \,l) \cr &= f_L (l) \prod\limits_{i = 1}^n \bigg[ {\displaystyle{{{\rm \beta}_i} \over {{\rm \eta}_i}} \left( {\displaystyle{{s_i - {\rm \gamma}_i} \over {{\rm \eta}_i}}} \right)^{\!{\rm \beta}_i - 1} \exp \bigg( { - \bigg( {\displaystyle{{s_i - {\rm \gamma}_i} \over {{\rm \eta}_i}}} \bigg)^{\!\!{\rm \beta}_i}} \bigg)} \bigg] .}$$

Then, according to Eq. (33), by integrating the joint PDF $f_{\bf Z} ({\bf z})$ in Eq. (35) in the failure region Ω defined in Eq. (34), system designers can obtain the probability of system failure.

In the second case, the components of ${\bf Z} = (S_1\comma \,S_2\comma \, \ldots\comma \,S_n\comma \,L)$ are normally distributed. From Eq. (31), the distribution of the reconstructed component limit-state function is $Y_i \sim N({\rm \mu} _{Y_i}\comma \, {\rm \sigma} _{Y_i} ^2 )$ , with ${\rm \mu} _{Y_i} = {\rm \mu} _{S_i} - {\rm \mu} _L $ , and ${\rm \sigma} _{Y_i} = \sqrt {{\rm \sigma} _{S_i} ^2 + {\rm \sigma} _L^2} $ , in which μ _L and σ _L are the mean and standard deviation of L. All the reconstructed limit-state functions ${\bf Y} = (Y_1\comma \,Y_2\comma \, \ldots\comma \,Y_n )$ then follow a multivariate normal distribution determined by the following mean vector and covariance matrix:

(36) $${\bf \mu} = \lpar {{\rm \mu}_{S_1} - {\rm \mu}_L\comma \, {\rm \mu}_{S_2} - {\rm \mu}_L\comma \, \ldots\comma \, {\rm \mu}_{S_n} - {\rm \mu}_L} \rpar \comma \,$$

(37) $$\sum = \left[ \matrix{ {\rm \sigma}_{Y_1}^2 \cr {\rm cov} (Y_2\comma \, Y_1 ) \cr \vdots \cr {\rm cov} (Y_n\comma \, Y_1 )} \quad \matrix{ {\rm cov} (Y_1\comma \, Y_2 ) \cr {\rm \sigma}_{Y_2}^2 \cr \vdots \cr {\rm cov} (Y_n\comma \, Y_2 )} \quad \matrix{ \cdots \cr \cr \ddots \cr \cdots} \quad \matrix{ {\rm cov} (Y_1\comma \, Y_n ) \cr {\rm cov} (Y_2\comma \, Y_n ) \cr \vdots \cr {\rm \sigma}_{Y_n}^2} \right] \comma \,$$

in which ${\mathop{\rm cov}} (Y_i\comma \, Y_j ) = {\rm \sigma} _L^2 $ . For this special case, because the distribution parameters of Y could be easily derived, it is more convenient to find p _fs using the PDF of Y, which is given by

(38) $$f_{\bf Y} ({\bf y}) = \displaystyle{1 \over {\sqrt {(2{\rm \pi} )^n \vert \sum \vert}}} \exp \left( { - \displaystyle{1 \over 2}({\bf y} - {\bf \mu} ) \sum\nolimits^{- 1} ({\bf y} - {\bf \mu} )} \right).$$

Then p _fs is easily obtained by integrating Eq. (38) in the failure region $\{ {\bf Y}\vert{\bf Y} \lt {\bf 0}\} $ .