An effective approach for probabilistic lifetime modelling based on the principle of maximum entropy with fractional moments

Accepted Manuscript

An effective approach for probabilistic lifetime modelling based on the principle of maximum entropy with fractional moments Xufang Zhang, Wei He, Yimin Zhang, Mahesh D. Pandey PII: DOI: Reference:

S0307-904X(17)30474-2 10.1016/j.apm.2017.07.036 APM 11886

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

30 October 2016 5 July 2017 18 July 2017

Please cite this article as: Xufang Zhang, Wei He, Yimin Zhang, Mahesh D. Pandey, An effective approach for probabilistic lifetime modelling based on the principle of maximum entropy with fractional moments, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.07.036

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • This paper presents a fractional moment method for probabilistic lifetime modelling of engineering systems. • Factional moments are calculated based on a small, simulated sample of remaining useful

CR IP T

life of the system. • The principle of maximum entropy (MaxEnt) with fractional moments is used to recover the system lifetime distribution.

• Applications of method are illustrated by several dynamical and discontinuous stochastic

AC

CE

PT

ED

M

AN US

systems.

1

ACCEPTED MANUSCRIPT

An effective approach for probabilistic lifetime modelling based on the principle of maximum entropy with fractional moments Xufang Zhanga,∗, Wei Hea,, Yimin Zhanga,, Mahesh D. Pandeyb, a

CR IP T

School of Mechanical Engineering and Automation, Northeastern University Shenyang, LN, China, 110819 b Department of Civil and Environmental Engineering, University of Waterloo Waterloo, ON, Canada, N2L 3G1

Abstract

The paper presents a fractional moment method for probabilistic lifetime modelling of uncer-

AN US

tain engineering systems. A novel feature of the method is the use of fractional moments, as opposed to integer moments commonly used in the structural reliability literature. The fractional moments are calculated from a small, simulated sample of remaining useful life of the system. And fractional exponents that are used to model the system lifetime distribution are

M

determined through the entropy maximization process, rather than assigned by an analyst in prior. Together with the theory of copula, the efficiency and accuracy of the proposed method

stochastic systems.

ED

are illustrated by the probabilistic lifetime modelling of several dynamical and discontinuous

Keywords: Fractional moments, Probabilistic lifetime modeling, System reliability analysis,

AC

CE

PT

Copula, Stochastic systems

∗

Corresponding Author; Tel/Fax:+86-24-83673819. Email address: [email protected] ( Xufang Zhang)

Preprint submitted to Applied Mathematical Modelling

July 25, 2017

ACCEPTED MANUSCRIPT

1. Introduction From the life-cycle perspective, probabilistic lifetime modelling of an engineering system provides important information for risk assessment of the system by evaluating the mean-time to failure, survival probability, dynamic hazard rate, and among others [1–4]. The key issue in the analysis is to fit a probability distribution to adequately model lifetime observations.

CR IP T

To achieve this, a prior distribution is judged empirically from available information, based on either on-site inspection data or simulation results of the system. Estimation of distribution parameters relies heavily on available mathematic techniques, such as the method of maximum likelihood estimation, statistical regression, and the method of moment. However, the bias and variation of the inference are highly sensitive to the assumed prior distribution [5].

AN US

An alternative approach to recover the lifetime distribution originates from the modern information theory. Jaynes [6] presented the MaxEnt as a rational approach to choose the most unbiased estimate amongst all possible candidates. The determined distribution function is consistent with inspection data and contains minimum spurious information. When only

M

moment constraints are specified, Shore and Johnson [7] proves that the MaxEnt is a uniquely correct method for the probabilistic inference that satisfies all consistency axioms. The readers

ED

are directed to a comprehensive review of the MaxEnt method in the literature [8]. Despite the conceptual elegance, practical difficulties of the MaxEnt approach began to

PT

emerge. The first issue is that a large number of integer moments (order > 4 ) are required to achieve a reasonable accuracy in modelling the distribution tail. If moments are analytically

CE

calculated, such as in random vibration analysis [9], it is not a major issue. However, the entropy maximization algorithm experiences numerical instability as the number of moment

AC

constraints becomes large [10]. Another thorny issue is that the tail of MaxEnt distribution becomes an oscillatory function due to non-monotonic nature of the polynomial embedded in the assumed entropy density function. Instead of a calculation procedure, if sample data are used to estimate higher moments for the MaxEnt analysis, then it poses an additional problem. The moments estimated form the data are known to have a large statistical error, which would creep into the MaxEnt procedure as well. Because of these practical limitations, the interest in engineering application of the MaxEnt procedure diminished over time. The difficulty in obtaining accurate estimate of lifetime distribution from small samples has been 3

ACCEPTED MANUSCRIPT

the main impediment to the application of MaxEnt in engineering risk assessment. Recently, the interest in MaxEnt has been highlighted with the emergence of fractional moments, i.e., real numbers (or fractions) of moment order instead of integers [11, 12]. The reason for the interest is that the distribution of a positive random variable can be characterized by using a finite number of fractional moments. In addition, a fractional moment

CR IP T

embodies statistical information about a large number of central moments as shown later in Section 2.1. An optimal set of fractional exponents are determined via the MaxEnt analysis instead of assigning prior values. Besides, an interesting feature of the fractional moment method is able to incorporate observational errors for inspection data modelling [13]. To sum up, shortcomings of the traditional MaxEnt method, as stated above, can be overcome by

AN US

using fractional moments as constraints in place of the regular integer moments.

Objective of the manuscript is to present a fractional moment method to quantify the uncertainty in lifetime prediction of engineering systems. The method is essentially an extension of conventional MaxEnt method based on integer moments. Together with the theory of copula, the advantage of using the fractional moment is illustrated by its effectiveness in

M

predicting lifetime distribution of stochastic and discontinuous engineering systems.

ED

The rest of the manuscript is organized as follows. Section 2 briefly introduces the motivation of using the fractional moment and its analytical results for commonly used lifetime randon variables. Section 3 introduces the use of fractional moments to recover lifetime dis-

PT

tribution of a system based on the principle of maximum entropy (MaxEnt). Together with remaining useful life prediction of a stochastic circuit model and lifetime analysis of a reactor

CE

core protection system, examples are presented in Section 4 to illustrate potential applications

AC

of the MaxEnt method in engineering practice. Conclusions are summarized in Section 5. 2. Fractional moments Considering a lifetime random variable T , an αth order fractional moment of T is mathe-

matically defined as α

E[T ] =

Z

tα f (t)d t

T

Here, α is a real constant, rather than an integer in conventional moments. 4

(1)

ACCEPTED MANUSCRIPT

Moments E[T α ] of the lifetime variable with negative fractions (i.e., α < 0) are referred to as the inverse moment problem in the literature [14]. A sufficient condition for the existence of fractional moments E[T α ] for α ∈ [−1, 1] is only if the underlying probability density function f (t) of the random variable is a bounded value when t → 0 [15]. For a positively defined and continuous lifetime random variable, some of its fractional moments are employed in the

CR IP T

following to recover the parent distribution f (t). 2.1. The need of fractional moments

Methods for recovering probability distribution of a lifetime random variable T from some of its integer moments have been extensively investigated in the literature [16]. The inverse

AN US

moment problem will become ill-conditioned when the large number of integer moments are employed for an accurate recovery. However, the instability of the integer moment problem can be overcome by compressing probabilistic information of the lifetime variable into the information carried by a few appropriately chosen, real-valued fractional exponents αk and the associated fractional moments E[X αk ]. Besides, it is easy to show that one order fractional

M

moment contains much more probabilistic information than its integer counterparts. This can be explained further by expressing a fractional moment as the weighted summation of many

ED

orders of integer moments in the following. Consider a fractional polynomial tα , its Taylor series expansion around a real constant c gives the result of

AC

CE

PT

∞   X α α−k t = c (t − c)k k k=0
in which, fractional binomial coefficients αk are defined as   α α(α − 1) · · · (α − k + 1) = k(k − 1) · · · 1 k α

(2)

(3)

wherein, the upper index α emphasizes that the coefficient makes sense when any real number
appears in the position[17]. For instance, one has −0.5 = (−0.5)(−1.5)(−2.5)/3! = −0.3125. 3 Assume that the lifetime random variable T has all orders of integer moments, i.e., moR ment integrals T |tk |f (t)dt are always bounded for k = 1, 2, · · · . By properly selecting the parameter c to guarantee the series convergence absolutely, the mathematical relationship 5

ACCEPTED MANUSCRIPT

between one order fractional moment and all central moments is derived as ∞   X  α α−k  α E[T ] = c E (T − c)k k k=0

(4)

By further expressing central moments of T by means of its integer moments: k

E (T − c)



=

k X

k−j

(−1)

  k k−j  j  c E T (j = 1, 2, · · · , k) j

(5)

CR IP T



j=0

It is clear to see that one order of the fractional moment can be expressed as the summation of all orders of conventional integer moments. Numerical results have shown that it requires a large number of integer moments (order > 50) to guarantee an accurate estimation of fractional

AN US

moments with the approximation [18, 19]. Based on the consideration, analytical results for the fractional moment E[T α ] associated with commonly used probability distributions are derived in the following.

2.2. Analytical results of fractional moments

M

We start to derive analytical results of the fractional moment based on the Lognormal distribution, and further extend to other commonly used lifetime distributions in later. The

PT

ED

probability density function (PDF) of a lognormally distributed random variable is defined as   1 (log t − ζ)2 f (t) = √ exp − (6) 2δ 2 2π δ t where δ (> 0) and ζ (> 0) represent scale and location parameters, respectively.

CE

Substituting for the expression, an αth order fractional moment defined in Eq.(1) deter-

AC

mines analytical result of the integration as   MTα = exp α(2ζ + αδ 2 )/2

(7)

Note that statistical moments of T have been explicitly determined as distribution parameters (δ and ζ) of the lifetime random variable and the fractional exponent α. Following the analytical result of fractional moment for the Lognormal distribution [20], Table 1 summarizes analytical results of the fractional moment for commonly used Weibull, Beta, Gamma and Chi-square distributions. 6

ACCEPTED MANUSCRIPT

Table 1: Analytical results of the fractional moment for commonly used lifetime random variables

PDF  2 f (t) = √2π1 tδ exp − (log2δt−ζ) 2
δ−1   f (t) = ζδ ζt exp − ( ζt )δ

Lognormal Weibull Beta Gamma

f (t) =

1 tζ−1 (1 B(ζ,δ)

f (t) =

1 tδ−1 ζ δ Γ(δ)

− t)δ−1

exp(− ζt )

δ

δ > 0, ζ > 0

Γ(δ)|Γ(α+ζ)| B(δ,ζ)|Γ(α+ζ+δ)|

δ > 0, ζ > 0

ζ α |Γ(α + δ)|/Γ(δ)

δ > 0, θ > 0

δ

f (t) = t 2 −1 exp(− 2t )/2 2 /Γ( 2δ )

Chi-square wherein Γ(x) =

R∞ 0

δ>0

exp(−z)z x−1 d z, and B(a, b) =

R1 0

4

Moment

Lognormal 3

2α |Γ(α + 2δ )|/Γ( 2δ )

AN US

4

ζ α |Γ(1 + αδ )|

z a−1 (1 − z)b−1 d z.

(b) 5

(a) 5

Weibull Beta

2

3

Gamma

2

Chi-Square

1

1 0 −1 −0.75 −0.5 −0.25

0

M

Moment

δ > 0, ζ > 0

Fractional Moment: MTα   exp α(2ζ + αδ 2 )/2

Parameter

CR IP T

Distribution

0.25 0.5 0.75

1

0 −1 −0.75 −0.5 −0.25

0

0.25 0.5 0.75

1

Fractional exponent α

ED

Fractional exponent α

Figure 1: Analytic fractional moments of various lifetime random variables (Fraction α ∈ [−1, 1] and distri-

PT

bution parameters δ = 0.5, ζ = 2.0)

CE

We assume that distribution parameters of the lifetime random variable T are δ = 0.5 and ζ = 2.0, respectively, and the fractional exponent α varies from −1 to 1. With results

AC

summarized in Table 1, analytical fractional moments are presented in Figure 1. It is noted that the fractional moment is a continuous function with respect to the real-

valued exponent α, rather than being limited to positive integers. If random samples are employed to calculate the fractional moment, an unbiased estimate of MTα will be: N 1 X α ˆ ti MT = N i=1

wherein, ti represents an ith realization of the lifetime random variable T . 7

(8)

ACCEPTED MANUSCRIPT

With the analytical result of fractional moments derived in the section, the principle of maximum entropy is further introduced to recover the true but unknown parent distribution of the lifetime random variable in the following. 3. Recovering lifetime distribution from some of its fractional moments

CR IP T

The inverse problem of deriving an unknown parent distribution from an array of integer moments has no unique solutions[21]. Researchers have resorted to heuristics a parametric distribution as a prior, and the assumption on the distribution type is a rather contentious one from both philosophical and practical points of view. The reasons are that tail probabilities tend to be highly sensitive to a parametric form, and assigning a parametric priori implies

AN US

adding spurious information to the inference.

An alternative approach for the objective originates from the modern information theory. Jaynes [6] presented the MaxEnt as a rational approach to choose the most unbiased estimate amongst all possible candidates. It is shown that the MaxEnt is a uniquely correct method for

M

the probabilistic inference that satisfies all consistency axioms. When only moment constraints are specified, the determined MaxEnt distribution is consistent with inspection data and

ED

contains minimum spurious information[7]. Contrary to conventional integer moment based entropy method in the literature, the paper proposes to use some of fractional moments to

PT

recover the lifetime distribution as shown in the following. 3.1. The MaxEnt lifetime distribution

CE

The entropy as one of probabilistic characteristic of a random quantity has been extensively

AC

documented in the literature: H[T ] = −

Z

T

  f (t) log f (t) dt

(9)

In the paper, we primary focus on parent distribution prediction of the lifetime random T

based on the principle of maximum entropy. Therefore, with specified constraints in terms of fractional moments of T , an estimate fˆ(t) of the true but unknown density function f (t) can

8

ACCEPTED MANUSCRIPT

be numerically determined by maximizing the Shannon entropy as follows:     Find: fˆ(t)    Z    ˆ Maximize: H[T ] = − fˆ(t) log fˆ(t) d t  T  Z     tαk fˆ(t)d t = MTαk (k = 1, · · · , m) Subject to:

(10)

T

k=0

in which, α0 = 0,

and λ0 = log

Z

exp

−

m X

(11)

  dt λk t αk

AN US

T



CR IP T

This leads to the MaxEnt density function fˆ(t) has a generic form of  X  m α k fˆ(t) = exp − λk t

T

T

k=1

Here, λ = [λ0 , · · · , λm ] and α = [α1 , · · · , αm ] represent Lagrange multipliers and fractional exponents of the MaxEnt distribution, respectively.

An interesting point in the formulation is that fractional exponents αk need not to be

M

specified a priori, rather they can be determined as a part of the entropy maximization process . Note that Eq.(10) is a typical moment constrained optimization problem. Alternatively,

ED

MaxEnt parameters can be optimized based on the Kullback-Leibler (K-L) divergence as follows [22]:

Z

   f (t) log f (t) fˆ(t) d t ZT Z     = f (t) log f (t) d t − log fˆ(t) f (t)d t

CE

PT

K[f, fˆ] =

T

(12)

T

Substituting for the MaxEnt density function fˆ(t) in Eq.(11), the K-L divergence can be

AC

further rewritten in a compact form as K[f, fˆ] = −H[T ] + λ0 +

m X

λk MTαk

(13)

k=1

wherein, the quantity H[T ] is the true entropy of the lifetime variable that is independent of parameters λ and α. Therefore, the minimization of the K-L divergence implies to minimize of the following function: I(λ, α) = K[f, fˆ] + H[T ] = λ0 + 9

m X k=1

λk MTαk

(14)

ACCEPTED MANUSCRIPT

To sum up, the unknown MaxEnt parameters α and λ can be determined by numerically implementing an unconstrained minimization problem in the following:    {λk , αk } (k = 1 · · · , m) Find: Z  X   X m m αk   Minimize: I(λ, α) = log exp − d t + λ t λk MTαk k  T

k=1

(15)

k=1

ment the unconstrained optimization problem.

CR IP T

A numerical routine based on the simplex search method [23] was programmed to impleWith the MaxEnt density function fˆ(t), numerical results for reliability and hazard rate

AN US

functions of the lifetime variable T can be easily derived as follows 
R  ˆ = ∞ exp − Pm λk τ αk dτ  Reliability (or survival) function: S(t)  k=0 t 
 P  αk exp − m λ k=0 k t ˆ
Hazard rate function: h(t) = R ∞ Pm αk dτ  t exp − k=0 λk τ   
 R∞ P  Cumulative hazard rate function: H(t) αk ˆ λ τ dτ = − log t exp − m k k=0

(16)

Note that the reliability function examines the probability that breakdowns of the system

M

occur beyond a given point in time, and the hazard rate function models the occurrence of only one, namely the first failure event, whereas H(t) indicates the accumulation of the hazard over

ED

the time period τ ∈ [0, t]. Therefore, to monitor the lifetime of a system across the support of its lifetime distribution, the hazard rate function is preferably used to summarize the survival

PT

probability based on inspection data [24].

CE

3.2. Example: Failure time data modelling of an electrical appliance To illustrate engineering applications of the MaxEnt procedure for probabilistic modelling

AC

of real engineering data, the example considers lifetime modelling of an electrical appliance based on 60 observations [25]. Note that a substantial number of small failure times implies that hazard rate function of the electrical appliance is relatively high for small realizations of the useful life.

The MaxEnt method with three-order fractional moments estimated based on 60 observations of the electronic appliance is employed to recover underlying lifetime distribution of the appliance, whereas benchmark result provided by the Nelson-Aalen estimator is used to verify the hazard rate function. Note that the Nelson-Aalen estimator does not require any 10

ACCEPTED MANUSCRIPT

distributional assumptions to estimate the hazard rate function [26]. Statistical regression test determines the Weibull distribution is the best known distribution to model the failure time data. Results for reliability function and hazard rate function of the electrical appliance are presented in Figure 2. −4 0

10

4

−2

0 0

2000 4000 6000 Failure time (cycles)

−3

M

ED PT

0.5

CE

Hazard function

ME−FM Weibull Data

1

0

2000 4000 6000 Operation time (cycles)

8000

(b) Reliability (Survival) Function

x 10

2000 4000 6000 Operation time (cycles)

5 4

(c) Hazard Rate Function

ME−FM Weibul Data

3 2 1 0 0

8000

AC

0 0

10

8000

(a) Probability Distribution

1.5

−1

10

AN US

2

ME−FM Weibull Data

CR IP T

6

Survival function

Data Weibull ME−FM

Cumulative hazard function

Density function

8

x 10

2000 4000 6000 Operation time (cycles)

8000

(d) Cumulative Hazard Function

Figure 2: Results for failure time data modelling of an electronic appliance

It is observed that results based on the Weibull distribution overestimate survival probabilities of the electronic appliance. Especially, it determines a monotone hazard function due to the subjective bias introduced by the predefined distribution type. The method of MaxEnt with three-order fractional moments, however, is able to capture the nonmonotone 11

ACCEPTED MANUSCRIPT

characteristic of the hazard rate function. The close agreement between MaxEnt results and observational data highlights the effectiveness of utilizing fractional moments for probabilistic modelling of engineering lifetime data. Similar results based on other failure dataset also validate effectiveness of the fractional moment method for lifetime data modelling, but are

CR IP T

not reported at here for the sake of brevity. 4. Probabilistic lifetime modelling of complex systems

Compared to probabilistic modelling of failure dataset in forgoing section, it is much more important to predict lifetime distribution of an uncertain system based on its physically mathematical model T = M(X) in terms of system useful life T and input random vector X.

AN US

Note that the vector X = [X1 , X2 , · · · , Xn ]T comprises all uncertain variables and the model function M(·) is seldom of explicitly available in engineering practices.

In this paper, random variables Xi (i = 1, · · · , n) are assumed to follow a known probability distribution. However, some of elements in X might exhibit statistical dependency. Based on

M

the consideration, the prediction of lifetime distribution based on the probabilistic model T = M(X) is firstly illustrated by means of a stochastic lithium-ion battery model. In the example,

ED

random variables are assumed as independent, but numerical solutions of the system useful life need to evaluate a three-order stochastic differential equation. The second example considers

PT

uncertain lifetime modelling of a reactor core protection system. Contrary to continuous performance function of the lithium-ion battery model, system lifetime assessment of the

CE

protection system requires to deal will a series min-max discontinuous functions. To account for various correlation structures between random variables, the copula theory is introduced

AC

for lifetime modelling of uncertain systems with dependent information. 4.1. Probabilistic remaining useful life prediction of a lithium-ion battery The examples considers remaining useful life prediction of a lithium-ion battery based on

a stochastic electricity model. The lithium-ion battery is used to power an unmanned aerial vehicle for test and prognostic [27]. The goal of the stochastic analysis is to predict the end-of-discharge time of the battery, which is an indicative measurement of remaining flight

12

ACCEPTED MANUSCRIPT

time of the unmanned aerial vehicle. Therefore, the end-of-discharge time of the battery is synonymous to the end-of-life of the system.

Rsp

Rs

ib

isp

Csp

is

Cs

Rp

V

CR IP T

Cb

I

ip

Figure 3: Equivalent circuit model of the lithium-ion battery[27]

AN US

An equivalent circuit model of the lithium-ion battery is depicted in Figure 3, wherein, a large capacitance Cb holds the charge capacity of the battery. Besides, its state-of-charge (SOC) is a time-dependent quantity[28]:

qmax − qb (t) Cmax

(17)

M

zSOC (t) = 1 −

in which, qb represents current charge in the battery, whereas qmax and Cmax denote its maximal

ED

possible charge and capacity, respectively.

The resistance related to surface overpotential is a nonlinear regression function with respect to zSOC :

PT

 Rsp (t) = r0 + r1 exp r2 [1 − zSOC (t)]

(18)

CE

as well as the effective capacitance Cb :

AC

2 3 Cb (t) = c0 + c1 zSOC (t) + c2 zSOC (t) + c3 zSOC (t)

(19)

Note that ri and ci are fitting parameters as summarized in Table 2.

Parameters Values

Table 2: Fitting parameters in Eqs.(18) and (19) [27]

r0

r1

0.0272

1.087 × 10−16

r2

c0

34.64 19.80

13

c1

c2

1745.00 −1.50

c3 −200.20

ACCEPTED MANUSCRIPT

Therefore, charges qk (t) (k = b, sp, s) of the battery model can be expressed as a third-order differential equation:

CR IP T

  (t)   + Cqsp + Cqss(t) −I q˙b (t) = − Cqbb(t)  Rp Rp sp Rp  
qs (t) qsp (t) 1 qb (t) 1 q ˙ (t) = − + − Cs R p + I sp C R C R R sp p sp b p   
  q˙s (t) = qb (t) − qsp − qs (t) 1 + 1 + I Cb R p Csp Rp Cs Rp Rs

(20)

Note that state variables of the circuit model are: the charge qb (t) in Cb , the charge qsp (t) in Csp , and the charge qs (t) in Cs . Therefore, substituting for numerical solutions of state variables, the time-variant output voltage V (t) of the lithium-ion battery can be obtained as qb (t) qsp (t) qs (t) − − Cb Csp Cs

AN US

V (t) =

(21)

It is of interest to predict the end-of-discharge(EOD) time of the battery, which is corresponding to the time of output voltage V (t) decreasing to a voltage threshold VEOD : (22)

M

TRUL (X) = arg {t : V (X; t) = VEOD }

Note that remaining useful life (RUL) TRUL (X) of the lithium-ion battery is an implicit

ED

function with respect to input vector X. With probabilistic characteristics of X in Table 3, the brutal Monte-Carlo simulation method based on a large number of samples is capable of

PT

providing empirical distribution of TRUL (X), which is referred to as benchmark results in the following.

CE

With initial values qb (0) = 31000 C, qsp (0) = qs (0) = 0 C and loading current I = 35 A, state variables of the circuit model and the RUL TRUL of the battery are deterministically

AC

calculated. Note that results in Figure 4 are obtained based on nominal mean-values of the input vector X. It is observed that output voltage V (t) of the circuit is a monotonically decreasing quantity with respect to discharging time of the lithium-ion battery. Given a realization of the voltage threshold VEOD , corresponding result of TRUL can be numerically obtained by numerical solving Eq.(21). For example, if VEOD = 10V, the RUL of the lithiumion battery would be 424 seconds as shown in Figure 4(b). Therefore, given a sample of random vector X, the numerical model allows to calculate corresponding sample of TRUL of the lithium-ion battery. 14

ACCEPTED MANUSCRIPT

Table 3: Probability distribution of random variables in the lithium-ion battery model

Distribution

Mean-value

COV

Rs

Ohm

Normal

0.005

0.05

Rp

Ohm

Normal

10000

0.05

Cs

Farad

Normal

115

0.05

Csp

Farad

Normal

316

0.05

Cmax

Farad

Normal

qmax

Coulomb

Normal

I

Ampere

Normal

qb (0)

Coulomb

Lognormal

42500

0.08

20100

0.08

35

0.10

31000

0.05

AN US

× 104

CR IP T

Unit

20

State variable: qb(t) 2

Time (s)

State variable: qsp(t)

M

200 0 40

Time (s) State variable: q (t) s

0

200

400

600

12

VEOD = 10V 8

ED

0

PT

20

Output voltage (V )

16

0 400

Charge (C)

Charge (C)

Charge (C)

4

Random variable

TRUL = 424s

4

800

Time (second)

0

200

400

600

800

Time (second)

(b) Remaining useful life (RUL)

CE

(a) Three state variables

AC

Figure 4: State variables and an illustration of RUL of the lithium-ion battery model

To predict overall probability distribution of the system RUL, an input vector of eight

random variables is simulated. Given this input, a sample of output voltage trajectory is calculated by numerically solving governing equation in Eqs.(20) and (21). A sample of RUL corresponding to the end-of-discharge criteria in Eq.(22) is obtained. Figure 5 presents empirical RUL distribution of the lithium-ion battery based on 105 samples. Note that the discharging voltage threshold is assumed as VEOD = 10V in the analysis. 15

ACCEPTED MANUSCRIPT

Table 4: Parameters for the MaxEnt density function of system RUL

Entropy

k

0

1

2

3

5.6807

λk

340.12

493.82

−337.42

14.008

αk

−−

−0.31904

0.09257

0.42579

CR IP T

Alternatively, fractional moments of TRUL estimated based on 500 samples are employed to predict parent distribution of the system RUL. Note that samples of X for the fractional moment calculation are generated based on the Helton sequence, since it is a low-discrepancy sequence compared with results simulated with random sampling methods [29]. With samples of the RUL, the procedure of MaxEnt with three-order fractional moments was employed

AN US

to estimate overall probability distribution of RUL. To evaluate accuracy of the proposed method, empirical distribution of TRUL estimated based on 105 are presented in Figure 5. And parameters of the MaxEnt distribution, namely, fractional exponents (α) and Lagrange multipliers (λ), are summarized in Table 4. These parameters were calculated based on 500 samples of remaining useful life of the lithium-ion battery. A close agreement between ME-FM

M

and MCS results in Figure 5 confirms the effectiveness of the fractional moment method for

6

ED

remaining useful life prediction of the dynamical stochastic system. × 10-3

100 MCS ME-FM

MCS ME-FM

AC

Failure probability

2

CE

PDF

4

3

10-1

PT

5

300

10-2

10-3

10-4

1

0 100

2.0875 10 −2

500

700

10-5 100

900

300

500

700

900

Discharging time (seconds )

Remianing useful life (seconds)

(a) Probability distribution

(b) Time-variant failure probability

Figure 5: Probability distribution of RUL and system failure probability of the lithium-ion battery (RUL: Remaining useful lifetime)

16

ACCEPTED MANUSCRIPT

With the determined complementary cumulative distribution function of TRUL presented in Figure 5(b), it is easy to assess failure probability of the battery model with respect to a specific operation time of the system. For example, failure probability of the system is about 2.0875 × 10−2 after 300 seconds operation of the unmanned aerial vehicle, which is fairly close

to its Monte-Carlo simulated result 1.995 × 10−2 based on 105 samples.

CR IP T

To investigate the effect of sample size N on the convergence rate of the entropy obtained by using the MaxEnt procedure, an experiment for the sample size N varying from 50 to 500 is carried out. Since analytical entropy of the system RUL distribution is no longer available, crude Monte-Carlo simulation with 105 samples is assumed to provide benchmark result of the entropy.

AN US

6

ME-FM MCS

5.6

M

Entropy

5.8

ED

5.4

100

PT

5.2 50

150

200

250

300

350

400

450

500

The number of samples (N)

CE

Figure 6: The effect of sample size N on the convergence rate of the entropy

Given results of the entropy presented in Figure 6, the MaxEnt procedure with 50 samples

AC

fails to provide reliable estimate of the entropy. But after increasing N from 50 to 500, it is observed the convergence of the entropy with respect to the sample size N as reference to the benchmark result 5.6873. Based on the observation, the MaxEnt procedure with 500 samples are used to prediction lifetime distribution of the battery model. Note that three-order fractional moments are used in forgoing experiments. Figure 7 further investigates the effect of moment number m on the convergence rate of the entropy. In the experiment, each fractional moment was estimated based on 500 Helton samples. It is 17

ACCEPTED MANUSCRIPT

7.5 ME-FM MCS

6.5

6

5.5 1

2

3

4

CR IP T

Entropy

7

5

6

AN US

The number of fractional moments (m)

Figure 7: The effect of the moment number m on the convergence rate of the entropy

seen that typically three-order fractional moments is found to be sufficient to achieve a good convergence of the entropy.

M

4.2. Lifetime analysis of a reactor core protection system with dependent random variables A reactor core protection system is designed to immediately terminate function of a nu-

ED

clear energy system in emergency. By breaking chain reaction, a cooling system is setup to remove decay heat from reaction core. All nuclear power plants have some forms of the pro-

PT

tection system. Control rod is designed as series of metal rods that are quickly inserted into the reactor core to absorb neutrons and rapidly terminate the nuclear reaction. Therefore,

CE

making use of the core protection system will assure that if equipment operates as designed, then all anticipated operational occurrences or postulated accidents will result in acceptable

AC

consequences.

Probabilistic lifetime analysis of the reactor core protection system depends heavily on the

assumption that the components operate independently within the system. Lessons learned from past accidents reveal that ignore such sources of dependence can result in nonconservative estimate of system failure probability. Statistical modelling of multivariate random variables is a difficult task due to the complexity in correlation structure and the curse of dimensionality of input random variables. The theory of copula simplifies the problem by separating the 18

ACCEPTED MANUSCRIPT

learning of marginal distributions from the learning of the multivariate correlation structure, and then linking them together into a dependent density model. RCPS000 +

RCPS002

+

•

RCPS003

CR IP T

RCPS001

A

RCPS004

RCF005-13

RPA300JA-RO

2/3

B

C

RPA300JA-RO

AN US

D

RCPS005

RCPS006

RCPS007

+

+

+

RCPS005-RC

RCPS006-MP

RCPS006-RC

RCPS007-MP

RCPS007-RC

E

F

G

H

I

J

ED

M

RCPS005-MP

PT

Figure 8: Result for fault-tree analysis of the reactor core protection system (Top event code: RCPS000)

The section considers probabilistic lifetime analysis of a reactor core protection system

CE

modelled by means dependent random variables. Figure 8 presents logical connection of sub-systems and components of the reactor core protection system. The fault-tree analysis

AC

determines the minimal cut-set of the system as  Minimum Cut Set := A, B, CD, EG, EH, F G, F H, EI, EJ, F I, F J, GI, GJ, HI, HJ (23)

wherein, {A, · · · , J} denote basic failure events of the core protection system.

Following the theory of reliability analysis, the system lifetime quantity TS (X) can be

19

ACCEPTED MANUSCRIPT

written as a logic function in terms of input random vector X:  TS (X) = min X1 , X2 , max{X3 , X4 }, max{X5 , X7 }, max{X5 , X8 }, max{X6 , X7 },

max{X6 , X8 }, max{X5 , X9 }, max{X5 , X10 }, max{X6 , X9 }, max{X6 , X10 }, (24) max{X7 , X9 }, max{X7 , X10 }, max{X8 , X9 }, max{X8 , X10 }

CR IP T

Note that the system lifetime function TS (X) is defined as a series minimum and maximal

connected functions, which cause conventional reliability methods, such as the first/second order reliability methods, are not applicable for system reliability analysis of the reactor core protection system [30, 31]. Moreover, reliability analysis of such a system always depends heavily on the assumption of independent variables within the system. However, it is more

AN US

realistic to consider statistical dependency among failure times of components. The paper proposes to use the theory of copula and the method of fractional moment to deal with reliability analysis of the reactor core protection system. The method is based on the concept of Monte-Carlo simulation but with a fairly small number of samples, which guarantees the

M

efficiency of the method for reliability analysis of complex systems. Table 5: Marginal lifetime distribution of random variables in the reactor core protection system

Event code

A

RCPS003

B

RCF005-013

C

Random variable

ED

Symbol

Distribution

Mean (×103 hours) COV

Weibull

8.2

0.15

X2

Lognormal

8.5

0.12

RPA300JA-RO

X3

Lognormal

9.8

0.22

D

RPA300JB-RO

X4

Weibul

10.5

0.18

E

RPS005-MP

X5

Lognormal

9.7

0.16

RPS005-RC

X6

Weibul

9.4

0.12

G

CE

PT

X1

RPS006-MP

X7

Lognormal

10.1

0.22

H

RPS006-RC

X8

Lognormal

9.9

0.21

I

RPS007-MP

X9

Weibul

10.6

0.18

J

RPS007-RC

X10

Weibul

8.9

0.19

AC

F

The copula model allows to learn marginal distributions separately from a multivariate dependence structure that links them together into a joint distribution function. It is one 20

ACCEPTED MANUSCRIPT

of most popular approaches for reliability modelling of a system with dependent random variables. Appendix A summarizes more details on the theory of copula. In the analysis, marginal distributions of random variables Xi are summarized as shown in Table 5. Moreover, to quantify the input correlation, lifetime probability distribution of components A and B are assumed as dependent variables (i.e., X1 and X2 ) based on the

CR IP T

Kendall’s correlation coefficient τ . Then, following the theory of copula, joint distribution function of random variables X1 and X2 are:

f12 (x1 , x2 ) = c(u1 , u2 ; θ) × f1 (x1 )f2 (x2 )

(25)

wherein, ui = Fi (xi ) and Fi (xi ) denotes marginal CDF of random variable Xi .

AN US

The Clayton, Frank and Gaussian copula models c(u1 , u2 ; θ) are summarized as shown in Table A.1, as well as the mathematical relationship between the Kendall’s correlation coefficient τ and the parameter θ for each copula model.

Table 6 summarizes numerical results of the copula parameter θ calculated based on vari-

M

ous combinations of the copula model c(u1 , u2 ; θ) and the Kendall’s correlation coefficient τ . Together with marginal PDFs and CDFs of Xi (i = 1, 2), Figure 9 depicts analytic results

ED

for the joint density function f12 (x1 , x2 ) and 500 samples of X1 and X2 simulated based on various copula models and correlation coefficients.

PT

Table 6: Realizations of parameter θ for various copula models and correlation coefficients

Kendall’s coefficient Calyton copula Frank copula 0.2222

0.9073

0.1564

τ = 0.5

2.0

5.7363

0.7071

τ = 0.9

18.0

38.2812

0.9877

AC

CE

τ = 0.1

Gaussian copula

It is observed that the dependent structure simulated by means of different copula models

are almost identical while T1 and T2 are weakly depended (i.e., τ = 0.1) as shown in Figure 9(a). However, for medium and strong statistical dependencies, the use of copula models is able to flexibly mimic various correlation configurations as presented in Figures 9(b) and 9(c). For example, the Calyton copula describes the correlation structure relatively tight at low quantile regions and loose for high quantiles. Most importantly, the Calyton model is able to 21

ACCEPTED MANUSCRIPT

14

Random variable -- T2

12 10 8 6 Data Calyton Copular

Data Frank Copular

2 2

4

6

8

10

12

14 2

Random variable -- T1

4

6

8

10

12

Random variable -- T1

Data Gaussian Copular

CR IP T

4

14 2

4

6

8

10

12

14

Random variable -- T1

(a) The Kendall’s correlation coefficient τ = 0.1 14

10 8 6 Data Calyton Copular

4

Data Frank Copular

2 4

6

8

10

12

Random variable -- T1

14 2

4

6

M

2

AN US

Random variable -- T2

12

8

10

12

14 2

Data Gaussian Copular

4

Random variable -- T1

6

8

10

12

14

Random variable -- T1

ED

(b) The Kendall’s correlation coefficient τ = 0.5

14

PT

Random variable -- T2

12 10

6 4

Data Calyton Copular

AC

2

CE

8

2

4

6

8

10

12

Random variable -- T1

Data Frank Copular

14 2

4

6

8

10

12

Data Gaussian Copular

14 2

4

Random variable -- T1

6

8

10

12

14

Random variable -- T1

(c) The Kendall’s correlation coefficient τ = 0.9

Figure 9: Joint distribution function of dependent variables X1 and X2 simulated by means of different copula models

22

ACCEPTED MANUSCRIPT

represent asymmetric joint density functions, whereas results for Frank and Gaussian copulas are almost symmetric. In general, the effect of a copula model on statistical configuration of joint density function f12 (x1 , x2 ) would be much more influential when a copula function combines with medium or strong levels of the statistical dependency, i.e, τ > 0.5. Based on the observation, the effect of

CR IP T

using different copula models on probabilistic lifetime prediction of the core protection system is further considered in the following. 4.2.1. The effect of copula models

To investigate the effect of using different copula models on probabilistic lifetime assess-

AN US

ment of the reactor core protection system, the Clayton, Frank and Gaussian models associated with a large correlation coefficient τ = 0.9 are considered.

Since analytical result for system lifetime distribution is no longer mathematically available, an experiment of Monte Carlo simulation based on various copula models has been implemented to estimate empirical entropy of system lifetime TS (X). To implement, a de-

M

pendent structure between X1 and X2 is modelled either by means of the Calyton, Frank, or Gaussian copula, whereas remaining random variables Xi (i = 3, · · · , 10) are assumed as

ED

independent. 5000 samples of input vector X are randomly generated, and corresponding realizations of system lifetime statistics TS (X) are obtained to calculate empirical probabil-

PT

ity distribution of TS (X). An estimate of its entropy, then, can be numerically obtained by substituting the empirical PDF into Eq.(9). After 200 rounds of the random experiment,

CE

simulation results of the entropy are presented as shown in Figure 10. It is observed that numerical estimates of the entropy based on various copula models are

AC

almost identical as presented in Figure 10(a), whereas the Gaussian copula is able to predict the entropy result with the least variation. In addition, the MaxEnt procedure based on 500 Helton sequences and three-order fractional moments has been implemented to predict the system lifetime distribution. Compared to results for the Monte-Carlo simulation, relative errors for the entropy estimated by means of the ME-FM method are presented in Figure 10(b). It is seen that the close agreement of entropy results with that of MCS further validates numerical accuracy of the MaxEnt method. 23

ACCEPTED MANUSCRIPT

2

1.43 ME-FM

1.4080

1.40

1.3988 1.3937

1.39 1.38 1.37

Calyton

Frank

0 -2 -4

CR IP T

Entropy

1.41

Relative error (× 10−3 )

1.42

-6 -8

Gaussian

Frank

Gaussian

(b) Relative errors of the entropy estimates

AN US

(a) Box-plot of entropy estimates

Calyton

Figure 10: Numerical results for the entropy of the system lifetime statistics TS (X) based on various copula models

M

10-1

ED

10-2

Benchmark Calyton Copula Frank Copula Gaussian Copula

10-3

CE

PT

Failure probability

100

2

4

6

8

10

12

Operation time (in thousand hours)

AC

10-4

Figure 11: Failure probabilities of the reactor core protection system predicted by means of the MaxEnt method based on various copula models

The cumulative lifetime distribution obtained by means of the MaxEnt method with threeorder fractional moment are pictured as shown in Figure 11. Three copula models (i.e., the Calyton, Frank and Gaussian copulas) were considered in the discussion to describe the statistical dependency between random variables X1 and X2 given the correlation coefficient 24

ACCEPTED MANUSCRIPT

τ = 0.9. It is interesting to see that the lifetime distribution predicted based on different copula models are almost identical, even though geometric configurations of f12 (x1 , x2 ) are different to some extents as demonstrated in Figure 9(c). Note that the system lifetime distribution curve is also useful to indicate the time-variant failure probability as reference to a specific operation time of the system. The effect of different correlation coefficients τ on the

CR IP T

system lifetime prediction is further investigated in the following. 4.2.2. The effect of correlation coefficients

The section investigates the effect of the Kendall’s correlation coefficient τ on time-variant reliability analysis of the reactor core protection system. To implement, the model of Gaussian

AN US

copula is assumed to describe statistical dependency between X1 and X2 . Following the random sample simulation procedure in Appendix B, the Monte-Carlo experiment based on the discontinuous system performance function in Eq.(24) is carried out. Figure 12 presents results of the entropy for various values of the Kendall’s correlation coefficient τ . It is seen that the entropy is positively correlated with the correlation coefficient

M

τ , which implies a bigger value of the Kendall’s correlation coefficient is, the larger uncertainty

ED

in system lifetime statistics TS (X) is. 1.5

AC

CE

Entropy

PT

1.4 1.3 1.2 1.1 1.0

-0.8 -0.6 -0.4 -0.2

0

0.2

0.4

0.6

0.8

The Kendall's correlation coefficient Figure 12: Box-plot of the entropy estimated based on various realizations of the Kendall’s correlation coefficient τ

25

ACCEPTED MANUSCRIPT

1.5 Data Regression

1.3 1.2 1.1

CR IP T

Entropy

1.4

y = 13.575x − 0.3349 R 2 = 0.9888

1.0 0.1

0.105

0.11

0.115

0.12

0.125

0.13

AN US

Coefficient of variation (COV)

Figure 13: Statistical regression model between COV and entropy of the system lifetime statistics (The Kendall’s correlation coefficient τ varies from −0.8 to 0.8 with an incremental step 0.2)

M

Figure 13 describes a linear regression model for a mathematical relationship between the entropy and coefficient of variation (COV) of TS (X) as the correlation coefficient τ varies from

ED

−0.8 to 0.8. Similar to results in Figure 13, the model further confirms the observation that the uncertainty and entropy of TS (X) are positively proportional to the Kendall’s correlation coefficient.

PT

Results for lifetime distribution of the reactor core protection system are numerically estimated by means of the MaxEnt with three-order fractional moments. As presented in

CE

Figure 14, PDFs of the system lifetime variable TS (X) modelled based on different correlation coefficients exhibit differences to some extents. But all reserve typically long LHS tails, which

AC

cause a throne accuracy issue of conventional reliability methods for reliability analysis of such critical systems. Figure 15 presents results for time-variant reliability of the reactor core protection sys-

tem estimated with the MaxEnt procedure based on 500 Helton sequences. It is interesting to observe that the system with a positive large correlation coefficient (e.g., τ = 0.9 in the example) exhibits a relatively higher system reliability level than systems with independent or negatively correlated random variables. This is because a positive correlation coefficient 26

ACCEPTED MANUSCRIPT

0.7 τ = −0.9 τ =0 τ = 0.9

0.6

0.4 0.3 0.2 0.1 2

4

6

8

10

12

AN US

0

CR IP T

PDF

0.5

System lifetime (in thousand hours) Figure 14: Lifetime distribution of the reactor core protection system estimated based on various correlation coefficients (Mean-values of TS (X) are 7.3896, 7.5286 and 7.7890 for τ = −0.9, 0 and 0.9, respectively)

ED

PT

System reliability

0.8 0.6

Benchmarks τ = -0.9 τ=0 τ = 0.9

M

1

CE

0.4

AC

0.2 0

3

5

7

9

11

13

Operation time (in thousand hours)

Figure 15: Time-variant reliability of the reactor core protection system modelled with various correlation coefficients

27

ACCEPTED MANUSCRIPT

can increase mean-value of TS (X) and accordingly safety margin of the system as observed in Figure 14. Note that a subjectively independent assumption (τ = 0) on input random variables might cause the nonconservative estimation on system failure probability, whereas a large negative valued correlation coefficient (e.g. τ = −0.9 in the example) is recommended for safety assessment of such a critical engineering system. In summary, the effect of cor-

CR IP T

relation coefficients τ has demonstrated much more signification implications than that of copula models for time-variant reliability analysis of the reactor core protection system with dependent random variables. 5. Conclusion

AN US

The paper presents an effective approach for probabilistic lifetime modelling of engineering systems based on the principle of maximum entropy (MaxEnt). The approach belongs to the method of moment in structural reliability theory, but using fractional moments instead of conventional integer moments so far in the literature. Advantages of the fractional moment

M

method is firstly illustrated by its ability of compressing probabilistic information of failure time observations into the information carried by a few real-valued fractional exponents,

ED

which are numerically optimized through the entropy maximization process. Engineering applications of the method are further demonstrated by probabilistic lifetime modelling of two

PT

complex systems: (a) the remaining useful life prediction based on a stochastic dynamical model, and (b) the time-variant reliability analysis of a reactor core protection system with

CE

dependent random variables. Note that analytical results for the probabilistic lifetime distribution of such systems are no longer available due to the implicit and multivariate system performance functions, the brutal Monte Carlo simulation based on 105 samples is assumed

AC

to provide benchmark results. In addition, the theory of copula is introduced to model joint distribution of dependent random variables. To improve sampling efficiency for the fractional moment calculation, the Helton sampling method is employed in the paper due to its low-discrepancy. Results have shown that typically three-order fractional moments (m = 3) estimated based on 500 Helton samples (N = 500) are found to be sufficient to achieve a good convergence of the entropy as reference to the MCS result. Besides, the system reliability predicted based on different copula models are almost identical, even though the 28

ACCEPTED MANUSCRIPT

corresponding joint density function are different to some extents. However, the correlation coefficient has exhibited significant influence on the results of system reliability analysis. Especially, the independent assumption of input random variables might cause nonconservative estimations of system failure probability, and a large negative valued correlation coefficient is recommended for risk assessment of critical engineering systems. Since the entropy method

CR IP T

uses only 500 Helton samples for the system reliability analysis, one can conclude that it is quite efficient compared to the brutal MCS method based on 105 samples. In summary, the fractional moment approach is generic and it provides an effective solution to probabilistic lifetime modelling of a wide class of engineering systems.

AN US

Acknowledgement

The authors would like to thank the subject editor and anonymous reviewers for their valuable comments and suggestions to improve the quality of the manuscript. XZ also appreciates Chinese National Natural Science Foundation (51405069), the Fundamental Research Funds for Central

financially supporting the research.

M

Universities (N150304009) and the Postdoctoral Science Foundation of China (2017M610182) for

ED

Appendix A. The theory of copula Let FX (x) denote the joint cumulative distribution function (CDF) of random vector X, and

PT

Fi (·) represents the marginal CDF of an ith random variable. A copula function C(u1 , · · · , un ) = Prob[U1 6 u1 , · · · , Un 6 un ] is a mapping C : [0, 1]n → [0, 1], wherein, Ui = Fi (ti ) are uniform

CE

random variables on [0, 1]. Then, the joint CDF FX (x) can be represented by the copula function

AC

C(u1 , · · · , un ) as

  FX (x1 , · · · , tn ) = C F1 (x1 ), · · · , Fn (xn )

(A.1)

Sklar [32] shows that for any joint distribution, such a copula function C(u1 , · · · , un ) satisfying

Eq.(A.1) always exists. Therefore, consider a bivariate vector X = [X1 , X2 ]T , the joint PDF of X can be defined by using the copula function as fX (x) = wherein, c(u1 , u2 ; θ) =

∂ 2 C(u1 , u2 ) ∂u1 ∂u2 = c(u1 , u2 ; θ)f1 (x1 )f2 (x2 ) ∂u1 ∂u2 ∂x1 ∂x2

∂ 2 C(u1 ,u2 ) ∂u1 ∂u2

(A.2)

is the copula density function modelled with parameter θ. And the

parameter θ can be specified in terms of the Kendall’s rank correlation parameter ρτ for most of

29

ACCEPTED MANUSCRIPT

bivariate copulas in the following: τ (θ) = 4

Z

0

1Z 1 0

C(u1 , u2 )c(u1 , u2 ; θ)du1 du2 − 1

(A.3)

Analytical results of the integration are summarized in Table A.1 for the Claytonm, Frank and Gaussian copulas. More details in regard to high-dimensional and other kinds of copulas, the readers

CR IP T

are directed to a comprehensive monograph in the literature [33]. Table A.1: Parameter θ for the commonly used copula functions

copula type

copula density function, c(u1 , u2 ; θ) the Kendall’s coefficient, τ (θ) 1+θ −θ (2+1/θ) (u1 u2 )1+θ (−1+u−θ 1 +u2 )

Clayton

Gaussian Independent

1

1 + 4[D(θ)/θ − 1]/θ

AN US

−θe−θ(u1 +u2 ) (e−θ −1) [(e−θ −1)+(e−θu1 −1)(e−θu2 −1)]2  ξ2 θ2 −2θξ1 ξ2 +ξ2 θ2  √ 1 exp − 1 2(1−θ2 ) 2 2 1−θ

Frank

θ/(θ + 2)

wherein ξi = Φ−1 (ui ) (i = 1, 2); and D(θ) =

2 arcsin(θ)/π 0

Rθ

t 0 et −1 dt.

parameter can be realized as τij =

M

Given n samples of random variables Xi and Xj , unbiased estimator for the Kendall’s correlation c−d c+d ,

wherein, c is the number of concordant pairs and d the

ED

number of discordant pairs in the sample.

PT

Appendix B. A Gaussian copula-based method for generating dependent samples Based on the Gaussian copula, the joint density function of random vector X = [X1 , · · · , Xn ]T

CE

can be defined as

fX (x1 , · · · , xn ) = c(u1 , · · · , un ; θ) ×

n Y

fi (xi )

(B.1)

i=1

AC

wherein, fi (ti ) is the marginal density function of an ith random variable, and c(u1 , · · · , un ; θ) is the n-dimensional Gaussian copula density function:

in which, ζ =



  c(u1 , · · · , un ; θ) = |θ|−1/2 exp − ζ T (θ −1 − I)ζ/2

(B.2)

T Φ−1 [F1 (x1 )], · · · , Φ−1 [Fn (xn )] is a uniform random vector, and Fi (xi ) is the

marginal cumulative distribution function of Xi , I is an n-dimensional identity matrix. Parameters θij ∈ [−1, 1] of the Gaussian copula density function are derived in Table A.1 as θij = sin(πρτ ij /2)

(for i, j = 1, · · · , n)

30

(B.3)

ACCEPTED MANUSCRIPT

wherein, ρτ ij are the Kendall’s correlation coefficient of random variables Xi and Xj (i, j = 1, · · · , n). To sum up, given the marginal CDF Fi (xi ) and the Kendall’s correlation matrix ρτ = [ρτ ij ]n×n , a simulation procedure for generating dependent samples based on the Gaussian copula can be summarized as follows: Step 1. Simulate a random vector x = [x1 , · · · , xn ]T from standard Normal distribution N (0, I);

CR IP T

Step 2. Evaluate θij = sin(πρτ ij /2), where ρτ ii = 1.0 and ρτ ij = ρτ ji for any i, j = 1, · · · , n; Step 3. Compute u = LT x, where L is the Cholesky decomposition of matrix θ = [θij ]n×n ;

Step 4. Obtain the sample of an ith lifetime random variable xi = Fi−1 [Φ(ui )], where Fi−1 (·) denotes the inverse CDF of random variable Ti if any i = 1, · · · , n;

AN US

Step 5. Repeat steps (1)-(4) for a large number of sample size N .

References

[1] J. Kim, W. Song, B. Kang, Stochastic approach to kinematic reliability of open-loop mechanism with

M

dimensional tolerance, Applied Mathematical Modelling 34 (5) (2010) 1225 – 1237. [2] H. Hong, N. C. Lind, Approximate reliability analysis using normal polynomial and simulation results, Structural Safety 18 (4) (1996) 329 – 339.

ED

¨ [3] B. C ¸ ekyay, S. Ozekici, Reliability, MTTF and steady-state availability analysis of systems with exponential lifetimes, Applied Mathematical Modelling 39 (1) (2015) 284 – 296.

PT

[4] H. Su, J. Li, Z. Wen, Z. Fu, Dynamic non-probabilistic reliability evaluation and service life prediction for arch dams considering time-varying effects, Applied Mathematical Modelling 40 (15-16) (2016) 6908 – 6923.

CE

[5] M. Pandey, Minimum cross-entropy method for extreme value estimation using peaks-over-threshold data, Structural Safety 23 (4) (2001) 345 – 363.

AC

[6] E. Jaynes, Information theory and statistical mechanics, Physical Review 108 (2) (1957) 171 – 190. [7] J. Shore, R. Johnson, Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy, IEEE Transactions on Information Theory 26 (1) (1980) 26 – 37.

[8] J. Kapur, H. Kesavan, Entropy Optimization Principles with Applications, Academic Press Inc., New York, 1992. [9] K. Sobezyk, J. Trebicki, Maximum entropy principle in stochastic dynamics, Probabilistic Engineering Mechanics 5 (3) (1990) 102 – 110.

[10] A. Tagliani, Hausdorff moment problem and maximum entropy: A unified approach, Applied Mathematics and Computation 105 (2-3) (1999) 291 – 305.

31

ACCEPTED MANUSCRIPT

[11] X. Zhang, M. D. Pandey, Structural reliability analysis based on the concepts of entropy, fractional moment and dimensional reduction method, Structural Safety 43 (3) (2013) 28 – 40. [12] M. D. Pandey, X. Zhang, System reliability analysis of the robotic manipulator with random joint clearances, Mechanism and Machine Theory 58 (12) (2012) 137 – 152. [13] E. E. Gomes-Gon¸calves, H. Gzyl, S. Mayoral, Density reconstructions with errors in the data, Entropy 16 (6) (2014) 3257 – 3272.

CR IP T

[14] A. Khuri, G. Casella, The existence of the first negative moment revisited, The American Statistician 56 (1) (2002) 44 – 47.

[15] W. Piegorsch, G. Casella, The existence of the first negative moment, The American Statistician 31 (9) (1985) 60 – 62.

[16] J. Ramberg, P. Tadikamalla, E. Dudewicz, E. Mykytka, A probability distribution and its uses in fitting data, Technometrics 21 (2) (1979) 201 – 214.

AN US

[17] R. L. Graham, D. E. Knuth, O. Patashnik, Concrete Mathematics: A Foundation for Computer Science, Addison-Wesley Publishing Company, Menlo Park, 1988, In: Chapter 5, Page 154, Equaiton 5.1. [18] P. Inverardi, A. Tagliani, Maximum entropy density estimation from fractional moments, Communications in Statistics-Theory and Methods 32 (2) (2003) 327–345.

[19] H. Gzyl, A. Tagliani, Hausdorff moment problem and fractional moments, Applied Mathematics and Computation 216 (11) (2010) 3319 – 3328.

M

[20] E. Taufer, S. Bose, A. Tagliani, Optimal predictive densities and fractional moments, Applied Stochastic Models in Business and Industry 25 (1) (2009) 57 – 71.

Arnold, London, 1994.

ED

[21] A. Stuart, J. Ord, Kendall’s Advanced Theory of Statistics: Distribution Theory, 6th Edition, Edward

[22] P. Inverardi, A. Tagliani, Maximum entropy density estimation from fractional moments, Communications

PT

in Statistics: Theory and Methods 32 (2) (2003) 327 – 345. [23] J. Lagarias, J. Reeds, M. Wright, P. Wright, Convergence properties of the Nelder-Mead simplex method

CE

in low dimensions, SIAM Journal on Optimization 9 (1) (1998) 112 – 147. [24] D. R. Cox, D. Oakes, Analysis of Survival Data, Chapman & Hal, London, 1984. [25] J. F. Lawless, Statistical Models and Methods for Lifetime Data, John Wiley & Sons Inc., New Jersey,

AC

2011.

[26] W. Nelson, Theory and applications of hazard plotting for censored failure data, Technometrics 14 (4) (1972) 945 – 966.

[27] S. Sankararaman, M. J. Daigle, K. Goebel, Uncertainty quantification in remaining useful life prediction using first-order reliability methods., IEEE Transactions on Reliability 63 (2) (2014) 603 – 619. [28] M. Chen, G. Rincon-Mora, Accurate electrical battery model capable of predicting runtime and i-v performance, IEEE Transactions on Energy Conversion 21 (2) (2006) 504 – 511. [29] L. Kocis, W. J. Whiten, Computational investigations of low-discrepancy sequences, ACM Transactions

32

ACCEPTED MANUSCRIPT

on Mathematical Software 23 (2) (1997) 266–294. [30] H. Madsen, S. Krenk, N. Lind, Methods of Structural Safety, Dover publications Mineola, New York, 2006. [31] G. Peri¸caro, S. Santos, A. Ribeiro, L. Matioli, HLRF-BFGS optimization algorithm for structural reliability, Applied Mathematical Modelling 39 (7) (2015) 2025 – 2035. [32] A. Sklar, Fonstions de r´epartition `a n dimensions et leurs marges, Publications de l’Institut de Statistique

CR IP T

de l’Universit´e de Paris 8 (1959) 229 – 231.

AC

CE

PT

ED

M

AN US

[33] R. B. Nelsen, An Introduction to Copulas, Springer-Verlag, Inc., New York, 2006.

33