# Time-Dependent Reliability-Based Design Optimization Utilizing Nonintrusive Polynomial Chaos

Authors: Yao Wang; Shengkui Zeng; Jianbin Guo
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Time-dependent reliability-based design optimization (RBDO) has been acknowledged as an advance optimization methodology since it accounts for time-varying stochastic nature of systems. This paper proposes a time-dependent RBDO method considering both of the time-dependent kinematic reliability and the time-dependent structural reliability as constrains. Polynomial chaos combined with the moving least squares (PCMLS) is presented as a nonintrusive time-dependent surrogate model to conduct uncertainty quantification. Wear is considered to be a critical failure that deteriorates the kinematic reliability and the structural reliability through the changing kinematics. According to  Archard’s wear law, a multidiscipline reliability model including the kinematics model and the structural finite element (FE) model is constructed to generate the stochastic processes of system responses. These disciplines are closely coupled and uncertainty impacts are cross-propagated to account for the correlationship between the wear process and loads. The new method is applied to an airborne retractable mechanism. The optimization goal is to minimize the mean and the variance of the total weight under both of the time-dependent and the time-independent reliability constraints.


---

## Body

## 1. Introduction

In recent decades, numerous endeavors have been made to develop the reliability-based design optimization (RBDO) methods, due to the fact that RBDO maximizes the performance under constrains of target reliability level accounting for various sources of uncertainty [1, 2]. These developed RBDO methodologies can be classified into two groups: time-independent RBDO and time-dependent RBDO.Time-independent RBDO, which assumes reliability constraints are time independent, has been studied widely [2]. However, in many engineering cases, because of the degradation and stochastic loads, systems deterioration with time is such a severe problem that it must be taken into account. For instance, wear is one of the most critical failures that substantially affect the life span of bearings, hinges, and other mechanisms, and it should be considered in design phase [3]. For these applications, time-dependent RBDO methods should be conducted because the reliability is time varying. Most literatures on this topic concentrate on the time-dependent structural reliability problem. Two basic scenarios are presented, which are out-crossing methods [4–7] and extreme value methods [8, 9].The out-crossing methods require the calculation of the crossing rate of the likelihood that the performance falls into the failure domain. The most fundamental equation of solving the out-crossing problems is the Rice formula. Kuschel and Rackwitz [4, 5] proposed a method based on the out-crossing approach to evaluate the time-dependent reliability in the context of first-order reliability methods (FORM) and asymptotic second-order reliability methods (SORM). The most probable point (MPP) in this method is defined as “point of maximum local crossing rate.” Rectangular wave renewal processes and Gaussian processes are derived for load models, and the optimization goal is to minimize the total cost which contains initial cost and failure cost. Streicher and Rackwitz [6] improved this method by considering the dependencies among different failure modes, and the numerical Laplace transforms are used for the treatment of aging components. Andrieu-Renaud et al. [7] presented an out-crossing-based method which is called PHI2. The limit state function is considered as a random process whose out-crossing of the zero level is to characterize.The extreme value methods consider the failure event to be equivalent to the event that the extreme value is greater than or less than the threshold in a time interval, and time-dependent reliability can be translated to a time-independent one if the distribution of the extreme value is identified. Li et al. [8] discovered that correlative information among the component random events is inherent in the equivalent extreme-value event. Thus an equivalent extreme value approach is presented to evaluate the structural reliability. In the approach proposed by Chen and Li [9], a virtual stochastic process associated to the extreme value of the studied stochastic process is constructed. The probability density evolution method is used to evaluate the instantaneous probability density function (PDF) of the virtual stochastic process, which would generate the PDF of the extreme value simultaneously.According to Bhatti [10], kinematic reliability is defined as the probability of output member’s position and orientation falling into a specified range from the desired position and orientation. Fewer researches focus on time-dependent kinematic reliability analysis. Zhang et al. [11, 12] proposed a mean value first-passage method based on the out-crossing methods in structural reliability to evaluate the time-dependent reliability of the function generator mechanisms. This method is under the assumption of normality for random dimension variables with small variances; thus the motion error is a nonstationary Gaussian process, but their work only contains the randomness of parameters. The degradation of mechanisms which is a significant factor to the time-dependent kinematic reliability was not taken into account.Wear is one of the most critical degradation failures for mechanisms. Worn joints would lead to motion errors as well as raise the stress. Thus wear would substantially deteriorate the kinematic reliability and structural reliability through the changing kinematics. Meanwhile, loads dominate the wear of joints, and worn joints would result in fluctuation of the motion of mechanisms and thus change the loads on joints conversely. The correlationship between the wear process and loads is required to be considered in the analysis of wear-related time-dependent reliability. Furthermore, each type of mechanisms has its particular wear characteristic which is determined by its function and operation condition. Because of the scarcity of the full-scale wear tests, it is hardly to obtain the accurate time-varying wear data in design phase especially for those newly designed mechanisms. Therefore, stochastic degradation process of system responses is difficult to obtain. All of these factors make the conventional out-crossing methods and extreme value methods inapplicable to the wear-related time-dependent RBDO problem.In this paper, a new time-dependent RBDO strategy is proposed to account for both of the time-dependent kinematic reliability and the time-dependent structural reliability as constrains. In order to perform the uncertainty quantification, a time-dependent surrogate model, which is called polynomial chaos combined with the moving least squares (PCMLS), is presented to approximate the stochastic process. Nonintrusive polynomial chaos (NIPC) is used to describe the stochastic nature of system responses at selected time points. Then moving least squares (MLS) method is employed to approximate the time variant functions of PCE coefficients. Because of the explicit polynomial formula and simple structure, Monte Carlo simulation (MCS) is able to be conducted on the surrogate model to evaluate the time-dependent kinematic reliability and structure reliability.To tackle with the scarcity of wear data, a multidiscipline reliability model is constructed according to Archard’s wear law. Disciplines in the model are closely coupled and uncertainty impacts are cross-propagated to simulate the wear process. The model takes an iterative simulation to generate the stochastic processes of system responses. The correlationship between the wear process and loads is considered in the model, and the uncertain inputs include random variables and stochastic loads. Thus the time-dependent degradation tests and data are not necessary.The proposed time-dependent RBDO method includes three stages: (1) multidiscipline reliability model and simulation for generation of degradation process; (2) PCMLS for time-dependent uncertainty quantification; (3) genetic algorithm (GA) for the global optimum.The whole approach is demonstrated at an airborne retractable mechanical system under the stochastic wind load. Wear of the hinge is considered to be the most critical failure. The multidiscipline reliability model, which involves kinematics model and structural FE model, is constructed. The optimization goal is to minimize the mean and variance of the total weight under the time-dependent probabilistic constraints related to the kinematic reliability and the structure reliability.The paper is organized as follows. In Section2, the theories of NIPC and MLS are presented, respectively, and the PCMLS model is derived. In Section 3, the theoretical basis and the numerical procedure of the multidiscipline reliability model are presented. In Section 4, the proposed time-dependent RBDO is formulated. The case of the airborne retractable mechanical system is studied in Section 5. Conclusions are drawn in Section 6.


## 2. Polynomial Chaos Combined with the Moving Least Squares

The polynomial chaos expansion is a promising surrogate model that uses a set of orthogonal polynomial bases to approximate the random space of the system response [13]. According to whether it requires the modification of the deterministic code, the polynomial chaos approaches can be divided into two groups: intrusive approach and nonintrusive approach. Intrusive approach calculates the unknown polynomial chaos coefficients by projecting resulting equations onto basis functions for different modes. It requires the modification of the deterministic code. Thus it is difficult, expensive, and time consuming for many complex computational problems [14]. On the contrary, in the nonintrusive PC (NIPC), simulations are used as black boxes and the calculation of chaos expansion coefficients is based on a set of simulation response evaluations. Most nonintrusive approaches are based on sampling or quadrature methods. Hosder et al. [14] applied the point-collocation NIPC to an aerospace problem with multiple uncertain variables. Sudret [15] proposed a nonintrusive regression-based PCE to conduct the global sensitivity analysis. Cheng and Sandu [16] proposed a least squares NIPC approach based on collocation at a low-discrepancy set of points which is demonstrated to have similar accuracy with the Galerkin approach by numerical experiments.Most NIPC methods are used to approximate the time-independent responses. However, a few literatures concentrate on the improvement of NIPC methods to propagate the time-dependent uncertainty. Witteveen et al. [17] proposed the probabilistic collocation for limit cycle oscillation (PCLCO) to modeling the long-term stochastic behavior of dynamical systems. PCLCO transforms the time-dependent issue into a time-independent one through the time-independent parametrization of the periodic response. Then the NIPC can be performed due to the independence of time, but PCLCO is only applicable to approximate the periodically time-dependent response. This paper proposes a surrogate model method which involves polynomial chaos combined with the moving least squares (PCMLS) to quantify more general time-dependent uncertainty. The surrogate model is essentially a time-dependent polynomial chaos expansion, and the coefficients of which are approximation functions achieved by the MLS method. The polynomial chaos expansion with time-dependent coefficients is used to approximate the stochastic process.


### 2.1. Polynomial Chaos Expansion

The polynomial chaos of a system response can be described as follows [18–20]:
(1)r(ξ)=a0+∑i1=1∞ai1Γ1(ξi1)+∑i1=1∞∑i2=1i1ai1i2Γ2(ξi1,ξi2)+∑i1=1∞∑i2=1i1∑i3=1i2ai1i2i3Γ3(ξi1,ξi2,ξi3)+⋯,
where r(ξ) is the random system response, ai is the coefficient of PCE, Γp(ξi1,…,ξip) is the polynomial of the selected base, ξi is a random variable with a suitable PDF agreed with the polynomial base, and p is the polynomial degree. According to the Wiener-Askey scheme, there are several families of polynomials (e.g., Hermite, Legendre, Jacobi, etc.) that could be employed depending on the random distributions. When ξi represents for a standard normal random variable, Hermite’s polynomials are selected as bases. The multidimensional Hermite polynomials in (1) is given by
(2)Γp(ξi1,…,ξip)=(-1)pe(1/2)ξTξ∂p∂ξi1⋯∂ξipe-(1/2)ξTξ,
where ξ is the vector of normal random variables {ξik}k=1p. In practical engineering, PCE contains limited input uncertainties. Thus (1) can be simplified as follows:
(3)r(ξ)=∑j=0Nc-1ajΨj(ξ),
where  Ψj(ξ)=∏i=1pψmij(ξi)=Γp(ξi1,…,ξip) and Nc is the total number of PCE coefficients, which can be calculated as
(4)Nc=1+n!(n-1)!+(n+1)!(n-1)!2!+⋯+(n-1+p)!(n-1)!p!=(n+p)!n!p!,
where n is the number of random variables in the system.The multidimensional Hermite polynomials form an orthogonal basis for the space of square-integrable PDFs, and the PCE is convergent in the mean-square sense [21]. In general, the approximation accuracy rises with the order of the PCE.Probabilistic collocation method is one of the most efficient NIPC methods. The coefficients of the PCE are obtained according to evaluations of the system response at selected collocation points, and these collocation points correspond to the roots of the polynomial of one degree higher than the order of the PCE [14]. Furthermore, the probabilistic collocation method just requires calls of the simulation the same as the number of PCE coefficients thus it is very efficient. The matrix form of probabilistic collocation method is describe as follows:
(5)[r(ξ1)r(ξ2)⋮r(ξNc)]=[Ψ0(ξ1)Ψ1(ξ1)⋯ΨNc-1(ξ1)Ψ0(ξ2)Ψ1(ξ2)⋯ΨNc-1(ξ2)⋮⋮⋱⋮Ψ0(ξNc)Ψ1(ξNc)⋯ΨNc-1(ξNc)]×[a0a1⋮aNc-1],
where ξj represents the jth set of collocation points. Total Nc sets of collocation points are required. Equation (5) can be written in a compact form
(6)r(Ξ)=Ψ(Ξ)a.
As a linear algebraic equation, (6) can be solved by Gauss elimination method and so forth. The sets of collation points are required to be chosen carefully to keep the matrix Ψ(Ξ) nonsingular, and the vector of coefficients can be derived as
(7)a=Ψ-1(Ξ)r(Ξ)=D(Ξ)r(Ξ).


### 2.2. Moving Least Squares

In order to reduce the cost of constructing numeral PCEs along the time period of degradation, the PCE is only built at some selected time points, and the moving least squares (MLS) method is employed to approximate the time-dependent functions of PCE coefficients. MLS is a generalization of the least squares technique, and it has become a widespread and powerful tool in interpolating and approximating implicit surfaces [22–25]. MLS starts with a weighted least squares formulation for an arbitrary fixed point in ℝd, and then this point moves over the entire parameter domain, where a weighted least squares fit is computed and evaluated for each point individually. Suppose a compact set Ω⊆ℝd is given and a continuous function f will be reconstructed from its values f(x1),f(x2),…,f(xN) on scattered pairwise distinct centers  X={x1,x2,…,xN}. Then the approximate value f^(x) is given by [23]
(8)minf^∈Πqd∑iw(∥x-xi∥)∥f^(xi)-f(xi)∥2,
where d represents the dimension of the function f, q represents the order of p(x) which is the polynomial basis vector, and w(∥x-xi∥) is the weighting function. The functional form and parameters of the weighting function determine the accuracy of the MLS approximation. The weighting function should satisfy the principles as follows [26].(1)
The weighting function is nonnegative,(9)w(∥x-xi∥)>0insidethesupportdomainΩ.(2)
The weighting function has compact support,(10)w(∥x-xi∥)>0outsidethesupportdomainΩ.(3)
The weighting function should be normalized,(11)∫Ωw(∥x-xi∥)dΩ=1.(4)
The weighting function should be monotone decreasing,   and the weight values should decrease with the increase of the distance fromx.The most common weighting functions are spline function, compactly supported radial basis function (CSRBF), Gaussian function, and so forth. Among these weighting functions, the spline functions are used widely because its order can be chosen to obtain high approximation accuracy. Note that the higher the order of the spline functions would not necessarily perform a better approximation, and the order of the spline function is decided by the highest order of the approximated function derivative.The approximation function can be written as(12)f^(x)=pT(x)b(x),x∈Ω,
where b(x) is the coefficient vector. For a one-dimension problem, the polynomial basis p(x) has the form
(13)pT=(1,x),q=1,pT=(1,x,x2),q=2.
It is defined as
(14)J=(pb-f)TW(x)(pb-f),
where f={f(x1),f(x2),…,f(xN)}T. In order to compute b(x), extreme value of J should be taken as
(15)∂J∂b=A(x)b(x)-B(x)f=0,
and b(x) can be derived as
(16)b(x)=A-1(x)B(x)f,
where
(17)A(x)=pTw(x)p,(18)B(x)=pTw(x).
Therefore, the approximation function can be derived as
(19)f^(x)=Φ(x)f,
where Φ(x)=pT(x)A-1(x)B(x).


### 2.3. The Proposed PCMLS Method

The multidiscipline reliability model is simulated step by step, and one step is defined as a task execution in which the mechanism would perform a required function. The PCE is constructed on system responses at some selected steps along the service life. After this, MLS method is employed to approximate the time-dependent functions of PCE coefficients. TotalNs tasks need to be completed which is decided by the length of service life and frequency of task arrival. Ti is defined as the ith selected step (or task), and m is the total number of selected steps that PCEs are required to be built on. According to (7), the coefficients for all of the PCEs are solved by the following equation:
(20)aT1=D(Ξ)rT1(Ξ),aT2=D(Ξ)rT2(Ξ),⋮aTm=D(Ξ)rTm(Ξ),
where rTi(Ξ) is the vector of responses at Ti and aTi is the vector of the coefficients of PCE at Ti. Equation (20) can be simplified as
(21)Λ=D(Ξ)R(Ξ),
where D is mNc×mNc matrix, R=[rT1,rT2,…,rTm]T is mNc×1 matrix, and Λ=[aT1,aT2,…,aTm]T is mNc×1 matrix.After the calculation ofΛ, a new matrix Η of PCEs coefficients is created. Each row represents coefficients of the PCE at a selected step, and each column represents a specific coefficient varying with steps. Η can be then expressed in a column form
(22)H=[a0,a1,…,aNc-1],
where a0 is the vector of the first coefficients of PCEs along with the steps. According to (19), the approximation functions of time-dependent PCE coefficients are derived as
(23)f^(T)=Φ(T)Η,
where f^(T)=[f^0(T),f^1(T),…,f^Nc-1(T)], f^i(T) is defined as approximation function of the ith coefficient, and T=[T1,T2,…,Tm] is the vector of selected steps. Then the conventional PCE can be translated into a time-dependent PCE, which is written as follows
(24)r(ξ,T)=∑j=0Nc-1f^j(T)Ψj(ξ).
The proposed PCMLS method is able to consider the statistical correlation between a pair of responses at two different time points. The multidiscipline reliability model is simulated step by step to describe the degradation behavior of the mechanical system. Uncertainties of the parameters and design variables are brought into the model at the beginning of the simulation. The system degradation at the current step is computed according to the system degradation at the step before. Thus these uncertainties would propagate through different disciplines as well as evolve over time, and the system responses at different time points have statistical correlationship. The statistical correlation information is contained in the samples of system responses at different time points.The proposed PCMLS method is a nonintrusive method to construct the time-dependent PCE. The PCMLS model uses samples of system responses to perform the NIPC at different time steps. Then the MLS is conducted to approximate the time-dependent PCE coefficients. In this manner, the statistical correlation information is extracted into the time-dependent coefficients of the PCE, which are actually functions of time, and these time-dependent PCE coefficients would make the polynomials of PCE vary with time and thus would make the statistical properties of the approximated responses vary with time as well. Therefore, the statistical correlation is considered in the PCMLS model.


## 2.1. Polynomial Chaos Expansion

The polynomial chaos of a system response can be described as follows [18–20]:
(1)r(ξ)=a0+∑i1=1∞ai1Γ1(ξi1)+∑i1=1∞∑i2=1i1ai1i2Γ2(ξi1,ξi2)+∑i1=1∞∑i2=1i1∑i3=1i2ai1i2i3Γ3(ξi1,ξi2,ξi3)+⋯,
where r(ξ) is the random system response, ai is the coefficient of PCE, Γp(ξi1,…,ξip) is the polynomial of the selected base, ξi is a random variable with a suitable PDF agreed with the polynomial base, and p is the polynomial degree. According to the Wiener-Askey scheme, there are several families of polynomials (e.g., Hermite, Legendre, Jacobi, etc.) that could be employed depending on the random distributions. When ξi represents for a standard normal random variable, Hermite’s polynomials are selected as bases. The multidimensional Hermite polynomials in (1) is given by
(2)Γp(ξi1,…,ξip)=(-1)pe(1/2)ξTξ∂p∂ξi1⋯∂ξipe-(1/2)ξTξ,
where ξ is the vector of normal random variables {ξik}k=1p. In practical engineering, PCE contains limited input uncertainties. Thus (1) can be simplified as follows:
(3)r(ξ)=∑j=0Nc-1ajΨj(ξ),
where  Ψj(ξ)=∏i=1pψmij(ξi)=Γp(ξi1,…,ξip) and Nc is the total number of PCE coefficients, which can be calculated as
(4)Nc=1+n!(n-1)!+(n+1)!(n-1)!2!+⋯+(n-1+p)!(n-1)!p!=(n+p)!n!p!,
where n is the number of random variables in the system.The multidimensional Hermite polynomials form an orthogonal basis for the space of square-integrable PDFs, and the PCE is convergent in the mean-square sense [21]. In general, the approximation accuracy rises with the order of the PCE.Probabilistic collocation method is one of the most efficient NIPC methods. The coefficients of the PCE are obtained according to evaluations of the system response at selected collocation points, and these collocation points correspond to the roots of the polynomial of one degree higher than the order of the PCE [14]. Furthermore, the probabilistic collocation method just requires calls of the simulation the same as the number of PCE coefficients thus it is very efficient. The matrix form of probabilistic collocation method is describe as follows:
(5)[r(ξ1)r(ξ2)⋮r(ξNc)]=[Ψ0(ξ1)Ψ1(ξ1)⋯ΨNc-1(ξ1)Ψ0(ξ2)Ψ1(ξ2)⋯ΨNc-1(ξ2)⋮⋮⋱⋮Ψ0(ξNc)Ψ1(ξNc)⋯ΨNc-1(ξNc)]×[a0a1⋮aNc-1],
where ξj represents the jth set of collocation points. Total Nc sets of collocation points are required. Equation (5) can be written in a compact form
(6)r(Ξ)=Ψ(Ξ)a.
As a linear algebraic equation, (6) can be solved by Gauss elimination method and so forth. The sets of collation points are required to be chosen carefully to keep the matrix Ψ(Ξ) nonsingular, and the vector of coefficients can be derived as
(7)a=Ψ-1(Ξ)r(Ξ)=D(Ξ)r(Ξ).


## 2.2. Moving Least Squares

In order to reduce the cost of constructing numeral PCEs along the time period of degradation, the PCE is only built at some selected time points, and the moving least squares (MLS) method is employed to approximate the time-dependent functions of PCE coefficients. MLS is a generalization of the least squares technique, and it has become a widespread and powerful tool in interpolating and approximating implicit surfaces [22–25]. MLS starts with a weighted least squares formulation for an arbitrary fixed point in ℝd, and then this point moves over the entire parameter domain, where a weighted least squares fit is computed and evaluated for each point individually. Suppose a compact set Ω⊆ℝd is given and a continuous function f will be reconstructed from its values f(x1),f(x2),…,f(xN) on scattered pairwise distinct centers  X={x1,x2,…,xN}. Then the approximate value f^(x) is given by [23]
(8)minf^∈Πqd∑iw(∥x-xi∥)∥f^(xi)-f(xi)∥2,
where d represents the dimension of the function f, q represents the order of p(x) which is the polynomial basis vector, and w(∥x-xi∥) is the weighting function. The functional form and parameters of the weighting function determine the accuracy of the MLS approximation. The weighting function should satisfy the principles as follows [26].(1)
The weighting function is nonnegative,(9)w(∥x-xi∥)>0insidethesupportdomainΩ.(2)
The weighting function has compact support,(10)w(∥x-xi∥)>0outsidethesupportdomainΩ.(3)
The weighting function should be normalized,(11)∫Ωw(∥x-xi∥)dΩ=1.(4)
The weighting function should be monotone decreasing,   and the weight values should decrease with the increase of the distance fromx.The most common weighting functions are spline function, compactly supported radial basis function (CSRBF), Gaussian function, and so forth. Among these weighting functions, the spline functions are used widely because its order can be chosen to obtain high approximation accuracy. Note that the higher the order of the spline functions would not necessarily perform a better approximation, and the order of the spline function is decided by the highest order of the approximated function derivative.The approximation function can be written as(12)f^(x)=pT(x)b(x),x∈Ω,
where b(x) is the coefficient vector. For a one-dimension problem, the polynomial basis p(x) has the form
(13)pT=(1,x),q=1,pT=(1,x,x2),q=2.
It is defined as
(14)J=(pb-f)TW(x)(pb-f),
where f={f(x1),f(x2),…,f(xN)}T. In order to compute b(x), extreme value of J should be taken as
(15)∂J∂b=A(x)b(x)-B(x)f=0,
and b(x) can be derived as
(16)b(x)=A-1(x)B(x)f,
where
(17)A(x)=pTw(x)p,(18)B(x)=pTw(x).
Therefore, the approximation function can be derived as
(19)f^(x)=Φ(x)f,
where Φ(x)=pT(x)A-1(x)B(x).


## 2.3. The Proposed PCMLS Method

The multidiscipline reliability model is simulated step by step, and one step is defined as a task execution in which the mechanism would perform a required function. The PCE is constructed on system responses at some selected steps along the service life. After this, MLS method is employed to approximate the time-dependent functions of PCE coefficients. TotalNs tasks need to be completed which is decided by the length of service life and frequency of task arrival. Ti is defined as the ith selected step (or task), and m is the total number of selected steps that PCEs are required to be built on. According to (7), the coefficients for all of the PCEs are solved by the following equation:
(20)aT1=D(Ξ)rT1(Ξ),aT2=D(Ξ)rT2(Ξ),⋮aTm=D(Ξ)rTm(Ξ),
where rTi(Ξ) is the vector of responses at Ti and aTi is the vector of the coefficients of PCE at Ti. Equation (20) can be simplified as
(21)Λ=D(Ξ)R(Ξ),
where D is mNc×mNc matrix, R=[rT1,rT2,…,rTm]T is mNc×1 matrix, and Λ=[aT1,aT2,…,aTm]T is mNc×1 matrix.After the calculation ofΛ, a new matrix Η of PCEs coefficients is created. Each row represents coefficients of the PCE at a selected step, and each column represents a specific coefficient varying with steps. Η can be then expressed in a column form
(22)H=[a0,a1,…,aNc-1],
where a0 is the vector of the first coefficients of PCEs along with the steps. According to (19), the approximation functions of time-dependent PCE coefficients are derived as
(23)f^(T)=Φ(T)Η,
where f^(T)=[f^0(T),f^1(T),…,f^Nc-1(T)], f^i(T) is defined as approximation function of the ith coefficient, and T=[T1,T2,…,Tm] is the vector of selected steps. Then the conventional PCE can be translated into a time-dependent PCE, which is written as follows
(24)r(ξ,T)=∑j=0Nc-1f^j(T)Ψj(ξ).
The proposed PCMLS method is able to consider the statistical correlation between a pair of responses at two different time points. The multidiscipline reliability model is simulated step by step to describe the degradation behavior of the mechanical system. Uncertainties of the parameters and design variables are brought into the model at the beginning of the simulation. The system degradation at the current step is computed according to the system degradation at the step before. Thus these uncertainties would propagate through different disciplines as well as evolve over time, and the system responses at different time points have statistical correlationship. The statistical correlation information is contained in the samples of system responses at different time points.The proposed PCMLS method is a nonintrusive method to construct the time-dependent PCE. The PCMLS model uses samples of system responses to perform the NIPC at different time steps. Then the MLS is conducted to approximate the time-dependent PCE coefficients. In this manner, the statistical correlation information is extracted into the time-dependent coefficients of the PCE, which are actually functions of time, and these time-dependent PCE coefficients would make the polynomials of PCE vary with time and thus would make the statistical properties of the approximated responses vary with time as well. Therefore, the statistical correlation is considered in the PCMLS model.


## 3. Multidiscipline Reliability Model Considering Wear

Wear is one of the most critical failures that substantially deteriorate the kinematic reliability and structural reliability of mechanisms in their service life. For a linkage mechanism, wear of the hinges would affect the sport stability, raise the stress in hinges, and lead to motion errors.Numerical wear prediction develops fast recently. It helps to study the effect of the wear evolution on the stress distribution and deformation of the motion pairs, which can hardly be measured with the experimental techniques. Söderberg and Andersson [27] presented a wear simulation method based on Archard’s wear law and adaptive FE method to simulate the wear of the brake pad under steady-state drag conditions. Rezaei et al. [3] proposed an adaptive wear modeling method to study the wear progress in radial sliding bearings. In this method, remeshing is performed both on the contact elements and their proximity elements. All of these proposed methods are demonstrated by experiments.However, these methods are based on the assumption that the sliding velocity and loads on the motion pairs hardly vary during the wear process. These two parameters are then considered to be boundary conditions, which are regarded as constants. In practice, the wear process would modify the kinematics continuously and lead to variation of the loads and the sliding velocity of motion pairs all the time. Traditional methods that neglect the variation of loads and the sliding velocity will cause inaccuracy at the wear simulating results. Meanwhile, in the above-mentioned methods, wear process is simulated at nominal state. Uncertainties of geometry and material are not taken into account.In this section, a multidiscipline reliability model that includes kinematics model, structural FE model, and wear model is constructed to simulate the wear process of the hinges in the linkage mechanism.


### 3.1. Wear Mechanism

To the plastically dominated wear, Archard’s law would serve as the appropriate model as discussed by Lim and Ashby [28]. In this model, the worn out volume is considered to be proportional to the normal load. The model is expressed mathematically as follows:
(25)dwdt=KFNvH,
where dw/dt  is defined to be the wear volume loss rate, K  is the dimensionless wear coefficient, FN is the applied normal force, H is the Brinell hardness of the softer material and it is usually considered to be a constant, and v denotes the relative sliding velocity. Since the wear depth is of interest, (25) is often written in the following form
(26)dhdt=kHPv,
where h is the wear depth, dh/dt  is the wear depth rate, and P is the contact pressure. We have kH=K/H,  and it is defined as the wear coefficient with the dimension of (Pa−1). The wear process is considered to be a time-dependent process.In practical engineering, the contact pressure varies with time because the real contact area changes during the wear process. Meanwhile, the contact pressure and relative sliding velocity vary with stochastic loadings. Thus both of these two parameters are time dependent. The numerical solution for the wear depth is obtained by estimating the differential form in (26) with a finite difference to yield the following updated formula:
(27)hi=hi-1+kHP(Θi,ti)v(Θi,ti)·Δt,
where hi is the wear depth at the ith step, Θi is the vector of stochastic loads at the ith step, and ti is ith step.


### 3.2. Multidiscipline Reliability Model

In Archard’s wear law, the wear rate is a function of both the pressure and the relative sliding velocity of two contact surfaces, and the pressure is obtained through nonlinear FEA, while the relative sliding velocity requires kinematics analysis. Meanwhile, the pressure and relative sliding velocity are varying with the wear process due to the fact that worn joints would change the kinematics. Therefore, the calculation of wear depth needs an iterative process.For a linkage mechanism, the multidiscipline reliability model includes kinematics model and structural FE model. The kinematics model is constructed on the entire mechanism to obtain the relative sliding velocity of hinges, and loads on hinges can also be given by the kinematics model. Then these loads are transmitted into the structural FE model, which is built with commercial finite element software ANSYS. FE model of hinges is solved by the nonlinear contact analysis to compute the stress distribution and contact pressure. After that, the conduct pressure and sliding velocity are delivered to the wear model to calculate the wear depth, according to which, the kinematics model and structural FE model would be updated in the next step, and the schematic of the multidiscipline reliability model is shown in Figure1.Figure 1
Schematic of the multidiscipline reliability model considering wear.In Figure1, Δhi, ΔΘi, Δvi, and ΔPi are the variation ranges of the wear depth, loads on hinges, relative sliding velocity, and contact pressure, respectively. Variation ranges of these parameters come from the dimension tolerances, parameter distributions, and stochastic environmental disturbance. Therefore, uncertainties exist in every discipline. As these disciplines are naturally closely coupled and uncertainty impacts are cross-propagated, the correlationship between the wear process and loads is considered in the multidiscipline reliability model.The simulation would generate the stochastic processes of stress and wear depths in hinges. The increase of wear depth would deteriorate the kinematic reliability, and the raise of stress in hinges would decrease the structural reliability. Thus the multidiscipline reliability model can be used to evaluate wear-related time-dependent kinematic reliability and structure reliability. The effect of wear process on the whole mechanism is taken into account by the coupling among wear model and other disciplines. The proposed method permits that values of parameters in wear model are obtained from relevant disciplines instead of determined by empirical assumptions; hence the multidiscipline reliability model is of higher accuracy. The simulation routine is described in Figure2.Figure 2
Simulation routine of the multidiscipline reliability model considering wear.In this procedure, there are two different types of discretizations which include (1) the continuous geometry discretized by finite elements; (2) the continuous material removal is approximated at a discrete set of steps. These discretizations are important factors to the accuracy of the simulation results. Thus the size of the discretizations is required to be determined carefully.Since the multidiscipline reliability model is built, wear process can be simulated and studied. Uncertainties exist in the model and propagate through different disciplines. MCS should be employed to quantify the time-dependent uncertainty. However, the non-linear contact FEA and kinematics analysis, which contain numerous nonlinear algebraic equations and differential equations required to solve numerically, are very computationally expensive. Furthermore, the iterative process in a single simulation leads to repeated calls of these expensive models. Therefore, traditional MCS are computational infeasible. The PCMLS model proposed in this paper can be employed to approximate the stochastic processes generated by the multidiscipline reliability model. Then the computational expensive implicit model is replaced by an explicit surrogate model with simple structure, and MCS is able to be conducted on the surrogate model which would substantially reduce the computational cost.


## 3.1. Wear Mechanism

To the plastically dominated wear, Archard’s law would serve as the appropriate model as discussed by Lim and Ashby [28]. In this model, the worn out volume is considered to be proportional to the normal load. The model is expressed mathematically as follows:
(25)dwdt=KFNvH,
where dw/dt  is defined to be the wear volume loss rate, K  is the dimensionless wear coefficient, FN is the applied normal force, H is the Brinell hardness of the softer material and it is usually considered to be a constant, and v denotes the relative sliding velocity. Since the wear depth is of interest, (25) is often written in the following form
(26)dhdt=kHPv,
where h is the wear depth, dh/dt  is the wear depth rate, and P is the contact pressure. We have kH=K/H,  and it is defined as the wear coefficient with the dimension of (Pa−1). The wear process is considered to be a time-dependent process.In practical engineering, the contact pressure varies with time because the real contact area changes during the wear process. Meanwhile, the contact pressure and relative sliding velocity vary with stochastic loadings. Thus both of these two parameters are time dependent. The numerical solution for the wear depth is obtained by estimating the differential form in (26) with a finite difference to yield the following updated formula:
(27)hi=hi-1+kHP(Θi,ti)v(Θi,ti)·Δt,
where hi is the wear depth at the ith step, Θi is the vector of stochastic loads at the ith step, and ti is ith step.


## 3.2. Multidiscipline Reliability Model

In Archard’s wear law, the wear rate is a function of both the pressure and the relative sliding velocity of two contact surfaces, and the pressure is obtained through nonlinear FEA, while the relative sliding velocity requires kinematics analysis. Meanwhile, the pressure and relative sliding velocity are varying with the wear process due to the fact that worn joints would change the kinematics. Therefore, the calculation of wear depth needs an iterative process.For a linkage mechanism, the multidiscipline reliability model includes kinematics model and structural FE model. The kinematics model is constructed on the entire mechanism to obtain the relative sliding velocity of hinges, and loads on hinges can also be given by the kinematics model. Then these loads are transmitted into the structural FE model, which is built with commercial finite element software ANSYS. FE model of hinges is solved by the nonlinear contact analysis to compute the stress distribution and contact pressure. After that, the conduct pressure and sliding velocity are delivered to the wear model to calculate the wear depth, according to which, the kinematics model and structural FE model would be updated in the next step, and the schematic of the multidiscipline reliability model is shown in Figure1.Figure 1
Schematic of the multidiscipline reliability model considering wear.In Figure1, Δhi, ΔΘi, Δvi, and ΔPi are the variation ranges of the wear depth, loads on hinges, relative sliding velocity, and contact pressure, respectively. Variation ranges of these parameters come from the dimension tolerances, parameter distributions, and stochastic environmental disturbance. Therefore, uncertainties exist in every discipline. As these disciplines are naturally closely coupled and uncertainty impacts are cross-propagated, the correlationship between the wear process and loads is considered in the multidiscipline reliability model.The simulation would generate the stochastic processes of stress and wear depths in hinges. The increase of wear depth would deteriorate the kinematic reliability, and the raise of stress in hinges would decrease the structural reliability. Thus the multidiscipline reliability model can be used to evaluate wear-related time-dependent kinematic reliability and structure reliability. The effect of wear process on the whole mechanism is taken into account by the coupling among wear model and other disciplines. The proposed method permits that values of parameters in wear model are obtained from relevant disciplines instead of determined by empirical assumptions; hence the multidiscipline reliability model is of higher accuracy. The simulation routine is described in Figure2.Figure 2
Simulation routine of the multidiscipline reliability model considering wear.In this procedure, there are two different types of discretizations which include (1) the continuous geometry discretized by finite elements; (2) the continuous material removal is approximated at a discrete set of steps. These discretizations are important factors to the accuracy of the simulation results. Thus the size of the discretizations is required to be determined carefully.Since the multidiscipline reliability model is built, wear process can be simulated and studied. Uncertainties exist in the model and propagate through different disciplines. MCS should be employed to quantify the time-dependent uncertainty. However, the non-linear contact FEA and kinematics analysis, which contain numerous nonlinear algebraic equations and differential equations required to solve numerically, are very computationally expensive. Furthermore, the iterative process in a single simulation leads to repeated calls of these expensive models. Therefore, traditional MCS are computational infeasible. The PCMLS model proposed in this paper can be employed to approximate the stochastic processes generated by the multidiscipline reliability model. Then the computational expensive implicit model is replaced by an explicit surrogate model with simple structure, and MCS is able to be conducted on the surrogate model which would substantially reduce the computational cost.


## 4. Time-Dependent RBDO Procedure

Time-dependent RBDO procedure is the organization of all the elements such as system optimization, system analysis, discipline analysis, and time-dependent uncertainty analysis. How to efficiently arrange these elements into an execution sequence is the key to realize time-dependent RBDO. This paper proposes a three-stage strategy, and the procedure is shown in Figure3.Figure 3
Procedure of time-dependent RBDO considering wear.In the first stage, multidiscipline reliability model considering wear is constructed. In the second stage, the PCMLS model is applied to approximate the stochastic processes. In the third stage, a double-loop optimization is employed. In the outer loop, optimization algorithm executes optimum search. At each iteration, the inter loop calls uncertainty analysis, which applies MCS to PCMLS model, to evaluate the design and its time-dependent uncertainty characteristics. Because the PCMLS model is an explicit formula both with high accuracy and simple structure, it would significantly improve the computational efficiency of the whole time-dependent RBDO procedure. The general time-dependent RBDO based on the PCMLS model is reformulated as follows:(28)minh~(x,d,q)s.t.P{(gi(T,ri(ξ,T),cfi)≤0)}≥RTDi,i=1,2,…,nTD,s.t.P{gj(x,d,q)≤0}≥RTIj,j=1,2,…,nTI,s.t.xL≤x≤xU,
where x is the vector of random design variables, d is the vector of deterministic parameters, and q is the vector of random parameters. ri(ξ,T) is the PCMLS model of the ith system response, cfi is the failure criterion of the ith system response, gi(T,ri(ξ,T),cfi) is the time-dependent limit state function, RTDi is the time-dependent reliability constrains, nTD is the total number of the time-dependent responses, gj(x,d,q) is the jth time-independent limit state function, RTIj is the time-independent reliability constrains, and nTI is the total number of the time-independent responses.


## 5. Case Study

The engineering case is an airborne retractable mechanical system. It is a four-link mechanism which consists of the hydraulic actuator, rods, and pin-and-lug hinges. It is designed to carry the functional device to move in accordance with the predetermined trajectory for switching between working and stopping positions.The system performs a retractable action when a task arrives, and it is required to work reliably during the service life. When the retractable system is performing, it is under the loads of wind and the weight of the functional device. The hinges are working in a nonlubricated environment. Therefore, wear of the hinges is considered to be a critical failure. Worn hinges give rise to fluctuation of the rod loads, decrease the kinematic accuracy, and lead to the mechanism lock during the movement. Furthermore, worn hinges would increase the stress in hinges which may cause them to fracture during the operation. The top hinge in Figure4 is the most dangerous one for it bears the largest load. Thus wear in this hinge is considered in the multidiscipline reliability model.Figure 4
Schematic of the airborne retractable mechanical system.


### 5.1. Multidiscipline Reliability Model

The multidiscipline reliability model is built to simulate the system degradation caused by wear according to the proposed method in Section3. The kinematics model of the mechanism and the FE model of the rods and hinges are built at the initial design point.Through the kinematics simulation, the kinematics parameters are calculated and shown in Figure5. The driving force is supplied by the hydraulic actuator to make the mechanism move properly. Because of the clearance between the pin and lug in hinges, there are fluctuations both in the driving force and instant relative velocity. The maximum driving force is 15423 N, which is below the rated pressure of the hydraulic actuator, and thus the mechanism would operate properly.Driving force and instant relative velocity between pin and lug computed by the kinematics model. (a) Driving force curve during the movement; (b) instant relative velocity curve during the movement.
(a)
(b)The FE model is built with commercial finite element software ANSYS. The materials of pin and lug are different, and the pin is more wear resistant than lug. Thus, the multidiscipline reliability model is based on the assumption that worn part only involves the lug. From Figure6(b), the contact mainly occurs at the head and shoulder of the lug where the wear may be more serious.FE model and nonlinear contact analysis results exported by ANSYS. (a) FE model of the top hinge; (b) contour picture of the pressure on the contact surface of lug.
(a)
(b)The process that the mechanism moves from the stopping position to the working position and then moves back is defined as a task. As the total time of each task is short, the wear loss of each task is very small. Therefore, each task is as a step, during which the kinematics analysis and FEA is conducted. At the end of the step, wear depth is calculated according to the parameters delivered by the kinematics model and FE model, and kinematics model and FE model would be updated in the next step according to the changing geometry. In the service life, the mechanism would act about 2000 times which is considered to be the total steps in the algorithm.The positions of three fixed supports and the length of the functional device arm are determined by requirements from higher level system. Thus they are considered to be boundary conditions in the optimization. The coordinate(XA,YA) defines the length and orientation of the upper rod and nether rod. Thus the coordinates of point A are significant design variables. Moreover, radiuses of the functional device arm and the upper rod contribute to the total weight of the mechanical system, and they are constrained by loads. Thus they are the other two design variables. In the multidiscipline reliability model, the parameter uncertainties involve machining errors and environment disturbance which are summarized in Table 1. The uncertain parameters of geometry are sampled at the beginning of the simulation, and the wind disturbance is sampled at each step.Table 1
Information on design variables and parameters.
Type
Symbol
Description
Mean
Std
Distribution
Random design variables
X
A
CoordinateX of point A
1176
2
Normal
Y
A
CoordinateY of point A
-
497
2
Normal
R
ur
Outer radius of the upper rod
30 mm
0.3 mm
Normal
R
arm
Outer radius of the device arm
50 mm
0.5 mm
Normal
Random parameters
L
sr
Length of actuator slider rod
1600 mm
3 mm
Normal
R
nr
Radius of the nether rod
20 mm
0.2 mm
Normal
G
d
Weight of device
2600 N
10 N
Normal
W
x
Wind disturbance ofX direction
400 N
5 N
Normal
W
y
Wind disturbance ofY direction
375 N
5 N
Normal
W
z
Wind disturbance ofZ direction
206 N
5 N
NormalAccording to results from the simulation at the initial design point, the time-dependent wear depth and contact pressure are shown in Figure7, while the time-dependent lug stress and maximum driving force are described in Figure 8.Time-dependent responses of wear depth and contact pressure. (a) Wear depth curve during the service life; (b) contact pressure curve during the service life.
(a)
(b)Time-dependent responses of stress of lug and maximum driving force. (a) Stress curve during the service life; (b) maximum driving force curve during the service life.
(a)
(b)The contact pressure in Figure7(b) is the average pressure of the whole contact surface. Thus the wear depth calculated by the contact pressure is also the average wear depth of lug. At the 2000th step, the wear depth is up to 13.4645 μm which is below the limit determined by the mechanism accuracy. With the increasing of the time step, the contact pressure shows a downward trend, and this trend makes the growth of the wear depth slows down gradually. These phenomena are consistent with the literatures [3, 21] which is demonstrated by experiments. The fluctuation of contact pressure is caused by the fluctuation of the driving force.Stress of lug raises as the wear depth grows, which illustrates that the wear would deteriorate the structure of lug. The clearance between pin and lug increases with the growth of the wear depth, which leads to the maximum driving force fluctuate with time. In every task, stochastic wind disturbance is loading on the device. Thus the stochastic wind disturbance is another reason to the vibration in stress and maximum driving force. Some experiments have been conducted on the prototype, and results are extracted to employ V&V on the multidiscipline reliability model.


### 5.2. The PCMLS Model

Figures7 and 8 illustrate that the wear depth and the stress of lug are time-dependent responses. The proposed PCMLS model is conducted to approximate these two responses.The preliminary sensitivity analysis (SA) is firstly employed to screen out the critical parameters whose deviation would have a great impact on the system responses. Total 10 parameters are considered in SA, which are shown in Table1. The forward difference method, which is commonly used in the preliminary sensitivity analysis, is executed due to its simplicity and rapidity. The sensitivity of the wear depth and von Mises stress of the lug are evaluated by several perturbation sizes of these 10 parameters to obtain reliable results.Based on the SA results and the practical engineering experiences, four parameters are chosen which include “coordinateX of point A,” “coordinate Y of point A,” “length of actuator slider rod,” and “weight of the device”. Thus the dimension of the PCE is four. Secondly, the sets of collocation points are selected, and the total number of sets is calculated by (4). Then collocation points of standard normal random variables are transformed into inputs of the multidiscipline reliability model. The transformation functions are listed in Table 2.Table 2
Transformation functions of critical random parameters.
Symbol
Description
Transformation functions
X
A
CoordinateX of point A
1176
+
2
ξ
1
Y
A
CoordinateY of point A
-
497
+
2
ξ
2
L
sr
Length of actuator slider rod
1600
+
3
ξ
3
G
d
Weight of the device
2600
+
10
ξ
4Thirdly, the sets of transformed values are brought into the multidiscipline reliability model to generate the time-varying system responses, respectively, and the responses at every 200 steps are selected to build the PCE, which means that the PCE are constructed at step 1, 200, 400, and so on till 2000. The coefficients of PCEs can be calculated by (20).After the construction of PCEs, MLS is employed to approximate the time-varying functions of PCE coefficients. The discrete points in Figure9 show that they have a nearly linear increase over time. In this case, as mentioned in Section 2.2, low-order spline weighting functions should be used to approximate the time-dependent functions of the PCE coefficients. In this case, the cubic spline function performed a better approximation. Thus the cubic spline function is chosen as the weighting function which is given as follows:
(29)w(x)={23-4x2+4x3x≤12,43-4x+4x2-43x312<x≤1,0x>1.
The order of PCMLS is determined by the PCE order. As 2-order PCE and 3-order PCE are performed to approximate the wear depth and stress of the lug at different time steps, MLS is used to obtain the time-dependent functions of the PCE coefficients, and the accuracy of PCMLS increases with the order. The comparative results are obtained at selected time steps which are 1, 500, 1000, 1500, and 2000 steps. The Monte Carlo simulation with 1000 samples is employed as a benchmark for the comparison. The statistics including the mean and the variance of the wear depth and stress of the lug are compared between the 2-order PCMLS and 3-order PCMLS, which are shown in Tables 3, 4, 5, and 6.Table 3
Comparative results of mean of the wear depth of the lug along the time axis.
Method
Samples
Mean of the wear depth of the lug (μm)
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.9413
5.5165
10.6144
13.2563
18.4850
2-order PCMLS
15
0.9450
5.5036
10.6283
13.1739
18.3849
Relative error
—
0.0039
0.0023
0.0013
0.0062
0.0054
3-order PCMLS
35
0.9388
5.5178
10.6205
13.2022
18.4134
Relative error
—
0.0027
0.0002
0.0006
0.0041
0.0039Table 4
Comparative results of the variance of the wear depth of the lug along the time axis.
Method
Samples
Variance of the wear depth of the lug
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.0042
0.0211
0.1501
0.1591
0.1843
2-order PCMLS
15
0.0052
0.0246
0.1349
0.1382
0.1625
Relative error
—
0.2381
0.1659
0.1013
0.1314
0.1183
3-order PCMLS
35
0.0045
0.0205
0.1409
0.1423
0.1677
Relative error
—
0.0714
0.0284
0.0613
0.1056
0.0901Table 5
Comparative results of the mean of the Von Mises stress of the lug along the time axis.
Method
Samples
Mean of the Von Mises stress of the lug (MPa)
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
40.7708
43.4974
46.6781
48.1653
50.4784
2-order PCMLS
15
40.6015
43.6492
46.4639
48.2551
50.3862
Relative error
—
0.0042
0.0035
0.0046
0.0019
0.0018
3-order PCMLS
35
40.6575
43.5669
46.6246
48.1785
50.5254
Relative error
—
0.0028
0.0016
0.0011
0.0003
0.0009Table 6
Comparative results of the variance of the Von Mises stress of the lug along the time axis.
Method
Samples
Variance of the Von Mises stress of the lug
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.4718
0.5598
0.5822
0.6020
0.6336
2-order PCMLS
15
0.3941
0.4821
0.4748
0.6624
0.6946
Relative error
—
0.1647
0.1388
0.1845
0.1003
0.0963
3-order PCMLS
35
0.4069
0.4889
0.6551
0.6593
0.6789
Relative error
—
0.1376
0.1267
0.1252
0.0952
0.0715Examples of MLS approximation curves for PCE coefficients of the wear depth and the stress of lug. (a) MLS approximation curve fora0 of the wear depth PCE; (b) MLS approximation curve for a0 of the stress of lug PCE.
(a)
(b)The accuracy of the PCMLS approximation is determined by two main factors: the PCE order and the number of MLS nodes. From Tables3, 4, 5, and 6, the PCMLS approximation accuracy rises with the PCE order. However, high-order PCE requires more samples, in other word, more execution of the simulation of the multidiscipline reliability model, and it would increase the computation burden intensively. Meanwhile, the MLS approximation accuracy, which rises with the number of nodes (in this case the samples of PCE coefficients at different time steps), would also have an influence to the accuracy of PCMLS. More MLS nodes would obtain high-accuracy approximation [26] but would introduce more computation at the same time. In order to apply the PCMLS to the real engineering problem, a tradeoff between accuracy and efficiency should be considered and determined.Figure9 shows two examples of MLS approximation curves for the PCE coefficients. It can be seen that the curves fit these points well, and the root-mean-square errors of these two approximation curves are 4.116e-7 and 0.3821, respectively. Approximation functions of the other 14 coefficients are also built by MLS. Then the PCMLS model is constructed completely, and the time consuming multidiscipline reliability model is replaced by the PCMLS model to conduct time-dependent uncertainty quantification using MCS.The wear depth of the hinge would deteriorate the kinematics accuracy of the mechanism. According to the requirements of the kinematics accuracy, the failure criterion of wear depth is determined as 50μm, and the safety strength of lug is 80 MPa. MCS is applied to the PCMLS model to calculate the time-dependent kinematic reliability and structure reliability with 10000 simulations. Note that the first passage time is a critical concept in the time-dependent reliability analysis, and we use the discretization of random processes and Monte Carlo simulation method to compute the time-dependent structure reliability. When the time-dependent limit state function first falls into the failure domain, this sample would be regarded as a failure and would be eliminated from the survival ones. The results are presented in Figure 10. The kinematic reliability after the 2000th task is 0.9890, and the structure reliability is 0.9850, which all satisfy the reliability constrains of the retractable mechanism.Time-dependent reliability computed through the PCMLS model. (a) Time-dependent kinematic reliability caused by the wear at the top hinge; (b) time-dependent structure reliability caused by the increase of stress of lug.
(a)
(b)


### 5.3. Time-Dependent RBDO

The time-dependent RBDO is formulated for the retractable mechanism. The design optimization problem is defined as to minimize the mean and the variance of the total weight subject to two time-dependent reliability constrains and three time-independent reliability constrains. The failures of retractable mechanism are defined as follows:(1)
the practical wear depth of lug exceeds the allowable wear depth determined by the kinematics accuracy;(2)
the stress in the lug is greater than the strength of the lug;(3)
the stress in the upper rod is greater than the strength of the upper rod;(4)
the stress in the device arm is greater than the strength of the device arm;(5)
the mechanism gets locked when the practical maximum driving force exceeds the rated output of the hydraulic actuator.The previous two failure criteria are related to the time-dependent reliability. The design problem is formulated as
(30)minF(μweight(x,d,q),σweight(x,d,q))s.t.P{g(T,rwd(ξ,T),cfwd)≤0}≥Rwd,P{g(T,rsl(ξ,T),cfsl)≤0}≥Rsl,P(gi(x,d,q)≤0)≥Ri,i=1,2,3,xL≤x≤xU,
where g(T,rwd(ξ,T),cfwd) is the time-dependent limit state function of the wear depth and g(T,rsl(ξ,T),cfsl) is the time-dependent limit state function of the stress of  lug. Also gi(x,d,q) is the time-independent limit state function of time-insensitive responses which include the maximum driving force of the mechanism, stress in the upper rod, and stress in the device arm.The RBDO in the case study is a multiobjective optimization problem which has mean and variance of the total weight as two objectives. The design objective is generally formulated as follows:(31)F=∑i=1l[(±)w1is1iμi+w2is2iσi],
where w1i and w2i are the weights and s1i and s2i are the scale factors for the mean and variance of the ith performance response, respectively.  l is the total number of the performance responses. The “+” sign is used when the response mean is to be minimized, and the “−” sign is be used when the response mean is to be maximized. This formulation belongs to the weighted sum method which is the most common approach to multiobjective optimization. The weights and scale factors should be positive to make the minimum of the objective to be Pareto optimal. In this case, we decide the scale factors to be the ideal points, which are shown as follows:
(32)s1=μweight∘=12,s2=σweight∘=0.1.
The weights of these two objectives are chosen as one, which means the mean and variance of the weight are equally important.Genetic algorithm (GA) is adopted in this paper because it is able to converge to the global solution rather than to a local solution. Since the two objectives are combined to form a single objective, a conventional single-objective GA is performed in the proposed RBDO. We have written the GA in MATLAB according to the literature [29] without modification.The optimal design results are shown in Table7, which involves the time-independent RBDO and the proposed time-dependent RBDO, and it is indicated that the optimal design results of these two methods are different. Because the proposed time-dependent RBDO considered the wear loss in the top hinge, its design result reallocated loads on three hinges of the mechanism, and this made the other two hinges partly share the load which was originally applied to the top hinge. Thus the wear in the top hinge was reduced to make the design of the mechanism satisfy the time-dependent reliability constrains. However, because of the increase of loads on the other two hinges, the upper rod and the device arm are required to bear more loads. Therefore, the outer radiuses of the upper rod and the device arm are greater than the design result of the time-independent RBDO. Thus the retractable mechanism in the design of the proposed time-dependent RBDO is a little heavier than the design of the time-independent RBDO.Table 7
Comparison of optimal design results between time-independent RBDO and proposed time-dependent RBDO.
Type
X
A
Y
A
R
ur (mm)
R
arm (mm)
Mean of total weight  (kg)
Variance of total weight
Time-independent RBDO
1235.01
-
421.40
22.45
45.24
15.3631
0.2383
Time-dependent RBDO
1239.28
-
415.93
22.85
47.19
15.9832
0.3369The time-dependent reliability curves by two methods are given in Figure11.Comparison of time-dependent reliability curves by time-dependent RBDO and time-independent RBDO. (a) Curves of kinematic reliability; (b) curves of structure reliability.
(a)
(b)Figure11 shows that the optimal design given by the proposed time-dependent RBDO is able to satisfy the time-dependent reliability constrains during the service life of 2000 tasks, but the design of time-independent RBDO would just be the optimal solution at the initial time. After the 1400th and 1800th tasks, the kinematic reliability and structural reliability fall below 0.98, respectively. The comparison indicates the necessity to account for the wear in hinges as well as time-dependent kinematic and structural reliability in design phase. The proposed time-dependent RBDO method is able to obtain the optimal design under the time-dependent reliability constrains. Thus it is applicable to engineering problem.


## 5.1. Multidiscipline Reliability Model

The multidiscipline reliability model is built to simulate the system degradation caused by wear according to the proposed method in Section3. The kinematics model of the mechanism and the FE model of the rods and hinges are built at the initial design point.Through the kinematics simulation, the kinematics parameters are calculated and shown in Figure5. The driving force is supplied by the hydraulic actuator to make the mechanism move properly. Because of the clearance between the pin and lug in hinges, there are fluctuations both in the driving force and instant relative velocity. The maximum driving force is 15423 N, which is below the rated pressure of the hydraulic actuator, and thus the mechanism would operate properly.Driving force and instant relative velocity between pin and lug computed by the kinematics model. (a) Driving force curve during the movement; (b) instant relative velocity curve during the movement.
(a)
(b)The FE model is built with commercial finite element software ANSYS. The materials of pin and lug are different, and the pin is more wear resistant than lug. Thus, the multidiscipline reliability model is based on the assumption that worn part only involves the lug. From Figure6(b), the contact mainly occurs at the head and shoulder of the lug where the wear may be more serious.FE model and nonlinear contact analysis results exported by ANSYS. (a) FE model of the top hinge; (b) contour picture of the pressure on the contact surface of lug.
(a)
(b)The process that the mechanism moves from the stopping position to the working position and then moves back is defined as a task. As the total time of each task is short, the wear loss of each task is very small. Therefore, each task is as a step, during which the kinematics analysis and FEA is conducted. At the end of the step, wear depth is calculated according to the parameters delivered by the kinematics model and FE model, and kinematics model and FE model would be updated in the next step according to the changing geometry. In the service life, the mechanism would act about 2000 times which is considered to be the total steps in the algorithm.The positions of three fixed supports and the length of the functional device arm are determined by requirements from higher level system. Thus they are considered to be boundary conditions in the optimization. The coordinate(XA,YA) defines the length and orientation of the upper rod and nether rod. Thus the coordinates of point A are significant design variables. Moreover, radiuses of the functional device arm and the upper rod contribute to the total weight of the mechanical system, and they are constrained by loads. Thus they are the other two design variables. In the multidiscipline reliability model, the parameter uncertainties involve machining errors and environment disturbance which are summarized in Table 1. The uncertain parameters of geometry are sampled at the beginning of the simulation, and the wind disturbance is sampled at each step.Table 1
Information on design variables and parameters.
Type
Symbol
Description
Mean
Std
Distribution
Random design variables
X
A
CoordinateX of point A
1176
2
Normal
Y
A
CoordinateY of point A
-
497
2
Normal
R
ur
Outer radius of the upper rod
30 mm
0.3 mm
Normal
R
arm
Outer radius of the device arm
50 mm
0.5 mm
Normal
Random parameters
L
sr
Length of actuator slider rod
1600 mm
3 mm
Normal
R
nr
Radius of the nether rod
20 mm
0.2 mm
Normal
G
d
Weight of device
2600 N
10 N
Normal
W
x
Wind disturbance ofX direction
400 N
5 N
Normal
W
y
Wind disturbance ofY direction
375 N
5 N
Normal
W
z
Wind disturbance ofZ direction
206 N
5 N
NormalAccording to results from the simulation at the initial design point, the time-dependent wear depth and contact pressure are shown in Figure7, while the time-dependent lug stress and maximum driving force are described in Figure 8.Time-dependent responses of wear depth and contact pressure. (a) Wear depth curve during the service life; (b) contact pressure curve during the service life.
(a)
(b)Time-dependent responses of stress of lug and maximum driving force. (a) Stress curve during the service life; (b) maximum driving force curve during the service life.
(a)
(b)The contact pressure in Figure7(b) is the average pressure of the whole contact surface. Thus the wear depth calculated by the contact pressure is also the average wear depth of lug. At the 2000th step, the wear depth is up to 13.4645 μm which is below the limit determined by the mechanism accuracy. With the increasing of the time step, the contact pressure shows a downward trend, and this trend makes the growth of the wear depth slows down gradually. These phenomena are consistent with the literatures [3, 21] which is demonstrated by experiments. The fluctuation of contact pressure is caused by the fluctuation of the driving force.Stress of lug raises as the wear depth grows, which illustrates that the wear would deteriorate the structure of lug. The clearance between pin and lug increases with the growth of the wear depth, which leads to the maximum driving force fluctuate with time. In every task, stochastic wind disturbance is loading on the device. Thus the stochastic wind disturbance is another reason to the vibration in stress and maximum driving force. Some experiments have been conducted on the prototype, and results are extracted to employ V&V on the multidiscipline reliability model.


## 5.2. The PCMLS Model

Figures7 and 8 illustrate that the wear depth and the stress of lug are time-dependent responses. The proposed PCMLS model is conducted to approximate these two responses.The preliminary sensitivity analysis (SA) is firstly employed to screen out the critical parameters whose deviation would have a great impact on the system responses. Total 10 parameters are considered in SA, which are shown in Table1. The forward difference method, which is commonly used in the preliminary sensitivity analysis, is executed due to its simplicity and rapidity. The sensitivity of the wear depth and von Mises stress of the lug are evaluated by several perturbation sizes of these 10 parameters to obtain reliable results.Based on the SA results and the practical engineering experiences, four parameters are chosen which include “coordinateX of point A,” “coordinate Y of point A,” “length of actuator slider rod,” and “weight of the device”. Thus the dimension of the PCE is four. Secondly, the sets of collocation points are selected, and the total number of sets is calculated by (4). Then collocation points of standard normal random variables are transformed into inputs of the multidiscipline reliability model. The transformation functions are listed in Table 2.Table 2
Transformation functions of critical random parameters.
Symbol
Description
Transformation functions
X
A
CoordinateX of point A
1176
+
2
ξ
1
Y
A
CoordinateY of point A
-
497
+
2
ξ
2
L
sr
Length of actuator slider rod
1600
+
3
ξ
3
G
d
Weight of the device
2600
+
10
ξ
4Thirdly, the sets of transformed values are brought into the multidiscipline reliability model to generate the time-varying system responses, respectively, and the responses at every 200 steps are selected to build the PCE, which means that the PCE are constructed at step 1, 200, 400, and so on till 2000. The coefficients of PCEs can be calculated by (20).After the construction of PCEs, MLS is employed to approximate the time-varying functions of PCE coefficients. The discrete points in Figure9 show that they have a nearly linear increase over time. In this case, as mentioned in Section 2.2, low-order spline weighting functions should be used to approximate the time-dependent functions of the PCE coefficients. In this case, the cubic spline function performed a better approximation. Thus the cubic spline function is chosen as the weighting function which is given as follows:
(29)w(x)={23-4x2+4x3x≤12,43-4x+4x2-43x312<x≤1,0x>1.
The order of PCMLS is determined by the PCE order. As 2-order PCE and 3-order PCE are performed to approximate the wear depth and stress of the lug at different time steps, MLS is used to obtain the time-dependent functions of the PCE coefficients, and the accuracy of PCMLS increases with the order. The comparative results are obtained at selected time steps which are 1, 500, 1000, 1500, and 2000 steps. The Monte Carlo simulation with 1000 samples is employed as a benchmark for the comparison. The statistics including the mean and the variance of the wear depth and stress of the lug are compared between the 2-order PCMLS and 3-order PCMLS, which are shown in Tables 3, 4, 5, and 6.Table 3
Comparative results of mean of the wear depth of the lug along the time axis.
Method
Samples
Mean of the wear depth of the lug (μm)
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.9413
5.5165
10.6144
13.2563
18.4850
2-order PCMLS
15
0.9450
5.5036
10.6283
13.1739
18.3849
Relative error
—
0.0039
0.0023
0.0013
0.0062
0.0054
3-order PCMLS
35
0.9388
5.5178
10.6205
13.2022
18.4134
Relative error
—
0.0027
0.0002
0.0006
0.0041
0.0039Table 4
Comparative results of the variance of the wear depth of the lug along the time axis.
Method
Samples
Variance of the wear depth of the lug
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.0042
0.0211
0.1501
0.1591
0.1843
2-order PCMLS
15
0.0052
0.0246
0.1349
0.1382
0.1625
Relative error
—
0.2381
0.1659
0.1013
0.1314
0.1183
3-order PCMLS
35
0.0045
0.0205
0.1409
0.1423
0.1677
Relative error
—
0.0714
0.0284
0.0613
0.1056
0.0901Table 5
Comparative results of the mean of the Von Mises stress of the lug along the time axis.
Method
Samples
Mean of the Von Mises stress of the lug (MPa)
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
40.7708
43.4974
46.6781
48.1653
50.4784
2-order PCMLS
15
40.6015
43.6492
46.4639
48.2551
50.3862
Relative error
—
0.0042
0.0035
0.0046
0.0019
0.0018
3-order PCMLS
35
40.6575
43.5669
46.6246
48.1785
50.5254
Relative error
—
0.0028
0.0016
0.0011
0.0003
0.0009Table 6
Comparative results of the variance of the Von Mises stress of the lug along the time axis.
Method
Samples
Variance of the Von Mises stress of the lug
1 step
500 step
1000 step
1500 step
2000 step
MCS
1000
0.4718
0.5598
0.5822
0.6020
0.6336
2-order PCMLS
15
0.3941
0.4821
0.4748
0.6624
0.6946
Relative error
—
0.1647
0.1388
0.1845
0.1003
0.0963
3-order PCMLS
35
0.4069
0.4889
0.6551
0.6593
0.6789
Relative error
—
0.1376
0.1267
0.1252
0.0952
0.0715Examples of MLS approximation curves for PCE coefficients of the wear depth and the stress of lug. (a) MLS approximation curve fora0 of the wear depth PCE; (b) MLS approximation curve for a0 of the stress of lug PCE.
(a)
(b)The accuracy of the PCMLS approximation is determined by two main factors: the PCE order and the number of MLS nodes. From Tables3, 4, 5, and 6, the PCMLS approximation accuracy rises with the PCE order. However, high-order PCE requires more samples, in other word, more execution of the simulation of the multidiscipline reliability model, and it would increase the computation burden intensively. Meanwhile, the MLS approximation accuracy, which rises with the number of nodes (in this case the samples of PCE coefficients at different time steps), would also have an influence to the accuracy of PCMLS. More MLS nodes would obtain high-accuracy approximation [26] but would introduce more computation at the same time. In order to apply the PCMLS to the real engineering problem, a tradeoff between accuracy and efficiency should be considered and determined.Figure9 shows two examples of MLS approximation curves for the PCE coefficients. It can be seen that the curves fit these points well, and the root-mean-square errors of these two approximation curves are 4.116e-7 and 0.3821, respectively. Approximation functions of the other 14 coefficients are also built by MLS. Then the PCMLS model is constructed completely, and the time consuming multidiscipline reliability model is replaced by the PCMLS model to conduct time-dependent uncertainty quantification using MCS.The wear depth of the hinge would deteriorate the kinematics accuracy of the mechanism. According to the requirements of the kinematics accuracy, the failure criterion of wear depth is determined as 50μm, and the safety strength of lug is 80 MPa. MCS is applied to the PCMLS model to calculate the time-dependent kinematic reliability and structure reliability with 10000 simulations. Note that the first passage time is a critical concept in the time-dependent reliability analysis, and we use the discretization of random processes and Monte Carlo simulation method to compute the time-dependent structure reliability. When the time-dependent limit state function first falls into the failure domain, this sample would be regarded as a failure and would be eliminated from the survival ones. The results are presented in Figure 10. The kinematic reliability after the 2000th task is 0.9890, and the structure reliability is 0.9850, which all satisfy the reliability constrains of the retractable mechanism.Time-dependent reliability computed through the PCMLS model. (a) Time-dependent kinematic reliability caused by the wear at the top hinge; (b) time-dependent structure reliability caused by the increase of stress of lug.
(a)
(b)


## 5.3. Time-Dependent RBDO

The time-dependent RBDO is formulated for the retractable mechanism. The design optimization problem is defined as to minimize the mean and the variance of the total weight subject to two time-dependent reliability constrains and three time-independent reliability constrains. The failures of retractable mechanism are defined as follows:(1)
the practical wear depth of lug exceeds the allowable wear depth determined by the kinematics accuracy;(2)
the stress in the lug is greater than the strength of the lug;(3)
the stress in the upper rod is greater than the strength of the upper rod;(4)
the stress in the device arm is greater than the strength of the device arm;(5)
the mechanism gets locked when the practical maximum driving force exceeds the rated output of the hydraulic actuator.The previous two failure criteria are related to the time-dependent reliability. The design problem is formulated as
(30)minF(μweight(x,d,q),σweight(x,d,q))s.t.P{g(T,rwd(ξ,T),cfwd)≤0}≥Rwd,P{g(T,rsl(ξ,T),cfsl)≤0}≥Rsl,P(gi(x,d,q)≤0)≥Ri,i=1,2,3,xL≤x≤xU,
where g(T,rwd(ξ,T),cfwd) is the time-dependent limit state function of the wear depth and g(T,rsl(ξ,T),cfsl) is the time-dependent limit state function of the stress of  lug. Also gi(x,d,q) is the time-independent limit state function of time-insensitive responses which include the maximum driving force of the mechanism, stress in the upper rod, and stress in the device arm.The RBDO in the case study is a multiobjective optimization problem which has mean and variance of the total weight as two objectives. The design objective is generally formulated as follows:(31)F=∑i=1l[(±)w1is1iμi+w2is2iσi],
where w1i and w2i are the weights and s1i and s2i are the scale factors for the mean and variance of the ith performance response, respectively.  l is the total number of the performance responses. The “+” sign is used when the response mean is to be minimized, and the “−” sign is be used when the response mean is to be maximized. This formulation belongs to the weighted sum method which is the most common approach to multiobjective optimization. The weights and scale factors should be positive to make the minimum of the objective to be Pareto optimal. In this case, we decide the scale factors to be the ideal points, which are shown as follows:
(32)s1=μweight∘=12,s2=σweight∘=0.1.
The weights of these two objectives are chosen as one, which means the mean and variance of the weight are equally important.Genetic algorithm (GA) is adopted in this paper because it is able to converge to the global solution rather than to a local solution. Since the two objectives are combined to form a single objective, a conventional single-objective GA is performed in the proposed RBDO. We have written the GA in MATLAB according to the literature [29] without modification.The optimal design results are shown in Table7, which involves the time-independent RBDO and the proposed time-dependent RBDO, and it is indicated that the optimal design results of these two methods are different. Because the proposed time-dependent RBDO considered the wear loss in the top hinge, its design result reallocated loads on three hinges of the mechanism, and this made the other two hinges partly share the load which was originally applied to the top hinge. Thus the wear in the top hinge was reduced to make the design of the mechanism satisfy the time-dependent reliability constrains. However, because of the increase of loads on the other two hinges, the upper rod and the device arm are required to bear more loads. Therefore, the outer radiuses of the upper rod and the device arm are greater than the design result of the time-independent RBDO. Thus the retractable mechanism in the design of the proposed time-dependent RBDO is a little heavier than the design of the time-independent RBDO.Table 7
Comparison of optimal design results between time-independent RBDO and proposed time-dependent RBDO.
Type
X
A
Y
A
R
ur (mm)
R
arm (mm)
Mean of total weight  (kg)
Variance of total weight
Time-independent RBDO
1235.01
-
421.40
22.45
45.24
15.3631
0.2383
Time-dependent RBDO
1239.28
-
415.93
22.85
47.19
15.9832
0.3369The time-dependent reliability curves by two methods are given in Figure11.Comparison of time-dependent reliability curves by time-dependent RBDO and time-independent RBDO. (a) Curves of kinematic reliability; (b) curves of structure reliability.
(a)
(b)Figure11 shows that the optimal design given by the proposed time-dependent RBDO is able to satisfy the time-dependent reliability constrains during the service life of 2000 tasks, but the design of time-independent RBDO would just be the optimal solution at the initial time. After the 1400th and 1800th tasks, the kinematic reliability and structural reliability fall below 0.98, respectively. The comparison indicates the necessity to account for the wear in hinges as well as time-dependent kinematic and structural reliability in design phase. The proposed time-dependent RBDO method is able to obtain the optimal design under the time-dependent reliability constrains. Thus it is applicable to engineering problem.


## 6. Conclusion

Existence of wear in mechanisms would deteriorate the time-dependent kinematic reliability and structural reliability through the changing kinematics. Thus wear-related time-dependent reliability should be taken into account in design phase. This work develops a time-dependent RBDO method considering wear as a critical failure. It involves both of the time-dependent kinematic reliability and time-dependent structural reliability as constrains.The PCMLS model that combines the NIPC with the MLS is presented to conduct the time- dependent uncertainty quantification. The NIPC is expanded to propagate time-dependent uncertainty with time-dependent coefficients achieved by MLS approximation functions.The multidiscipline reliability model, which includes kinematics model and structural FE model, is constructed according to Archard’s wear law to generate the stochastic processes of system responses. As these disciplines are closely coupled and uncertainty impacts are cross-propagated, the correlationship between the wear process and loads is considered in the model. The PCMLS model is applied to approximate the stochastic processes generated by the multidiscipline reliability model. Then MCS is conducted on the PCMLS model to evaluate the time-dependent kinematic reliability and structural reliability.The procedure of the three-stage time-dependent RBDO is given, and the new method is demonstrated at an airborne retractable mechanism. The optimization goal is to minimize the mean and the variance of the total weight. Constraints include both of the time-dependent and time-independent reliabilities. The optimal design is compared to the time-independent RBDO result. The comparison indicates that it is necessary to account for wear and time-dependent reliability in design phase, and the proposed time-dependent RBDO method is applicable to engineering problem.



---

*Source: 513261-2013-06-18.xml*