Modeling and analysis for time redundant systems with a given mission window

Accepted Manuscript Modeling and Analysis for Time Redundant Systems with a Given Mission Window He Yi, Lirong Cui, Jingyuan Shen PII: DOI: Reference:

S0360-8352(18)30516-3 https://doi.org/10.1016/j.cie.2018.10.036 CAIE 5475

To appear in:

Computers & Industrial Engineering

Received Date: Revised Date: Accepted Date:

28 March 2018 29 August 2018 16 October 2018

Please cite this article as: Yi, H., Cui, L., Shen, J., Modeling and Analysis for Time Redundant Systems with a Given Mission Window, Computers & Industrial Engineering (2018), doi: https://doi.org/10.1016/j.cie.2018.10.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.

Modeling and Analysis for Time Redundant Systems with a Given Mission Window He Yia, Lirong Cuib,* and Jingyuan Shenb a School b School

of Mathematics and Statistics, Beijing Institute of Technology, Beijing, China

of Management and Economics, Beijing Institute of Technology, Beijing China

Abstract Redundant systems take an important role in the field of reliability. Most of related researches have considered component redundancy, however, time redundancy which is also very common in practice hasn’t been discussed enough. In this paper, we assume that the time needed for a mission cannot be negligible and the mission must be completed in a given mission window which is long enough for several mission repeated tries. Based on this assumption, the main model is investigated together with its three special cases by the Laplace transform method, and their properties and relationships are discussed briefly. An extended model with different failure modes is taken into consideration by using the Monte Carlo method, and some extended discussions are given to compare the Laplace transform method with the Monte Carlo simulation method. The optimal analysis for cost and revenue are also presented for the main model. Then numerical examples are shown to illustrate results obtained in this paper. Finally, conclusions and discussions are presented.

Keywords: Time redundant system, mission window, mission repeated tries, Laplace transform method, Monte Carlo method

Nomenclature

[0,τ ]

given mission window.

τ

time that mission repeated tries are forced to stop.

Xi

ith mission repeated duration.

f ( x)

*

probability density function of X i , i = 1, 2, … .

Corresponding author. E-mail address: [email protected] (L. R. Cui). 1

F ( x)

distribution function of X i , i = 1, 2, … .

b

common upper bound of X i , i = 1, 2, … .

N (t )

number of mission repeated tries in interval [0, t ] for the main model, and

N (t ) = N j (t ) for Model j , j = 1, 2,3 . p

probability that the first long mission repeated try succeeds.

q

probability that the first long mission repeated try fails, and q = 1 − p .

a

threshold value such that mission repeated tries less than a must fail.

Yi

ith repair duration.

g (t )

probability density function of Yi , i = 1, 2, … .

G (t )

distribution function of Yi , i = 1, 2, … .

qki

failure probability of the ith mission repeated try in the main model if it is a long mission repeated try (whose duration is greater than a ) and qki = r i q . k

ki

number of long mission repeated tries before the ith mission repeated try.

r

parameter in the main model that describes the decreases in failure probability

qki after each failure of long mission repeated tries, and 0 < r < 1 .

H (t )

probability that the mission succeeds before time t for the main model and

H (t ) = H j (t ) for Model j , j = 1, 2,3 . h(t ), h j (t )

derivative functions of H (t ), H j (t ) with respect to t .

L(t ), L j (t )

partial derivative functions of H (t ), H j (t ) with respect to parameter q .

T

time to the mission success for the main model, and

H (τ )

mission success probability for the main model, and H (τ ) = H j (τ ) for Model j .

m

T = T j for Model j .

number of failure modes in the extended model.

2

X% i

ith mission repeat duration in the extended model and X% i = min{Z1i ,

Z 2i , … , Z mi } = Z ei ,i , (ei = 1, 2, … , m) .

Z vi

ith mission repeat duration under the vth (v = 1, 2,… , m) failure mode, whose probability density function is

f%v (t ) and distribution function is

F%v (t ) .

ei

failure mode of the ith mission repeated try.

Y%i

ith repair duration in the extended model, whose probability density function is g% v (t ) and distribution function is G% v (t ) if X% i = Z vi .

q%ei (i )

failure probability of the ith mission repeated try in the extended model if it k%

is a long mission repeated try and q%ei ( i ) = r%ei ei ( i ) q%ei , 0 < r%ei ≤ 1 .

k%v (i )

number of long mission repeated tries before the ith mission repeated try that have failure mode v (v = 1, 2,… , m) .

q%v

failure probability of the first long mission repeated try that has failure mode

v (v = 1, 2,… , m) .

r%v

parameter in the extended model that describes the decreases in failure probability after each long mission repeated try with failure mode

(v = 1, 2,… , m) , and 0 < r%v ≤ 1 . C (τ )

total cost in mission interval [0,τ ] .

R(τ )

net revenue in mission interval [0,τ ] .

f 0 , f1 , f 2

functions of H (τ ), EN (τ ) and q, r in the formula of ER (τ ) .

c

total funds for this mission in interval [0,τ ] .

3

v

B (c, d ; α , β )

Beta distribution with probability density function

xα −1 (d − x) β −1 I{c ≤ x ≤ d } d

∫u

α −1

(d − u )

β −1

.

du

c

1.

Introduction

Reliability of different kinds of systems has been extensively investigated in the literature (Ushakov, 1994; Barlow & Proschan, 1996; Lin et al., 2016; Zhu et al., 2017; Lin & Nguyen, 2018). For systems such as spacecraft, guaranteeing their high reliability is an extremely important issue, thus redundancy should be taken into consideration. As mentioned in Zhao et al. (2015), the most popular way of redundancy is to add redundant components to the system. If the redundant components are used to replace the failed component during the system operation, the system is called standby/cold redundant system; If the redundant components are used in parallel with the original components, the system is called active/hot redundant system (Eryilmaz, 2017; Zhao et al., 2017). Many systems familiar to us are active redundant systems, for example, parallel system in Barlow & Proschan (1996), k-out-of-n: F system in Linton & Saw (1974), consecutive-k-out-of-n: F system in Chang et al. (2000), mconsecutive-k, l-out-of-n system in Eryilmaz and Boushaba (2012), Cui et al. (2015) and so on. Warm standby redundancy, which balances the high operation cost of hot standby and the long restoration delay of cold standby, was also considered in researches (Rao & Gupta, 2000; Levitin et al., 2013; Jia et al., 2017). There are also other ways of redundancy, for example, functional redundancy, loadsharing redundancy and time redundancy. As introduced in Xu & Liao (2016), functional redundant systems are composed of multifunctional components which can be made to provide multiple functions. For example, different kinds of sensors are used to monitor different environmental parameters such as temperature, humidity, light intensity, pressure, wind speed, wind direction and so on. If more than one kind of these sensors are designed on one platform (for example, a silicon die), multiple environmental parameters will be monitored simultaneously, which means the platform can be regarded as a multifunctional component (Pallaro et al., 2004; Ivanov, 2012; Roozeboom et al., 2013; Nagel et al., 2013). Xu & Liao (2016) studied the reliability of a one-shot system containing multifunctional 4

components and formulated a reliability allocation problem with the objective of maximizing system reliability. Traditional component redundant systems waste resources by having all units operating at their highest performance levels, then load-sharing redundant systems, whose loads can be shared by all the functioning units, are introduced; The failure of a component will result in a higher load for the remaining components which therefore induces a higher failure rate in each of them (Shao & Lamberson, 1991; Liu, 1998; Kvam & Pena, 2005; Balakrishnan, 2011). Shao & Lamberson (1991) analyzed reliability and availability of a shared-load repairable k-out-of-n: G system under Markovian assumption by using inverse Laplace transform. Liu (1998) developed a proportional hazard model with Weibull baseline reliability to calculate reliability of a load-sharing k-out-of-n: G system with non-i.i.d. components. Inferences about parameters were also investigated for load-sharing systems in Kvam & Pena (2005), Balakrishnan (2011). Time redundant systems, which have been considered in Heidtmann (1982), Lisnianski et al. (2000), are systems whose total mission is a sequence of several phases and must be executed during constrained time. Denote the time needed to perform the system total mission as T and the constrained time as τ , and then reliability of the time redundant system can be given by the probability P{T ≤ τ } that the random variable T

is no more than the

constant/random variable τ (Lisnianski & Jeager, 2000). Lisnianski & Jeager (2000) further considered component redundancy for servers of a time redundant system under the assumption that failures of the system only lead to delay of mission performance rather than the destruction of the incomplete results of the system operation. As mentioned in Langer (2001), time redundancy can be unsupplemented (the time resource is determined before mission and cannot be changed afterwards), instantly supplemented (the time resource is determined after each repair) or gradually supplemented (the time resource increases during up periods, usually linearly). Wu & Hillston (2015) introduced two types of time redundant models with their mission reliability considered for semi-Markov systems. The first model required that the system should remain operational continuously for a minimum time within the given mission time interval and the second one 5

required that the total operational time of the system within the mission time window must be greater than a given value. In practice, there are many time redundant systems (Lisnianski et al., 2000). For example, control systems have a mission to drive some object from one point to another point and the mission should be completed in constrained time; communication systems have a mission to transmit information/material from one place to another and the transmission should be constrained in a given window. However, the exiting models are far from enough to describe so many kinds of time redundant systems in practice. In this paper, we investigate a new time redundant model together with its three special cases, and we also provide some extended discussions of the new model. Our main model can be applied to the spaceflight telemetry, tracking, and control (TT&C) system in Wu & Hillston (2015). Different from research like Wu & Hillston (2015), our models are not for performance of the TT&C system itself, but for its mission of data transmission. For a TT&C system, the data transmission from the spacecraft to the observation station on the ground must be accomplished in a given mission time window, during which the spacecraft is passing overhead the observation station when it orbits the earth. As shown in Figure 1, the solid circle means the earth and the dotted circle means the orbit of the spacecraft. Data transmission can only be done when the spacecraft orbits form point 1 to point 2 and the corresponding time is the mission window for data transmission of the TT&C system. In that time window, the data transmission task can be executed several times until it finally succeeds or until the time window ends, which means a time redundant system model will be needed to fit this practical example. (Figure 1 insert here about) For the TT&C system, speed of data transmission can be affected by many factors, for example, environment condition, network speed, performances of devices and so on. If the data transmission task is executed several times, obviously the time needed for each time, which we called mission repeated duration, may be different. According to the historical data about mission repeated durations, we can assume that they are independent identically distributed and then fit them with a proper distribution. However, the estimated distribution we find may have a larger value range than the true one. For this problem, we can estimate 6

the shortest time needed for a mission repeated duration by assuming that it happens at the highest transmission speed. If durations of some mission repeated tries are less than the shortest time we need, we can obviously deduce that they were interrupted by terrible influence of the environment (for example, sudden weather change or enemy countermeasures). Therefore, we call them short mission repeated tries and believe that they fail definitely. Long mission repeated tries can also fail, and after their failures, we can reduce the failure rate to some extent by finding out reasons and acting accordingly to avoid similar failures. Repair of devices can be carried out after each failure, no matter long ones or short ones. Time of repair can also be fitted by a proper distribution according to historical data. To guarantee high mission success probability in the given mission window, the time needed to complete a mission should be much less than the given mission window, which allows several mission repeated tries in this window. In our main model, we assume that the time durations of all mission repeated tries, which we call mission repeated durations, are independent and identically distributed random variables. Mission repeated tries fail definitely if their durations are less than or equal to a threshold value a , and mission repeated tries fail with decreasing probability if their durations are greater than a . We also consider that repair of device after each mission failure will cause delay of the next mission repeated try. By setting some parameters to be constant, there are three special cases of the main model discussed in this paper, denoted as Model 1, Model 2 and Model 3, respectively. To compare the Laplace transform method with the Monte Carlo method for numerical calculations of the main model, some extended discussions are also presented in this paper, including an extended model with different failure modes. The specific kind of models in this paper has not been investigated before, and can be applied to many different kinds of practical systems, which makes it a useful theoretical tool to investigate time redundant systems. The using of Laplace transform makes calculations for these models much more efficient than traditional methods such as the Monte Carlo method. Discussions on parameters are also presented to provide guidance on system design, and cost and revenue are also taken into consideration in some optimal problems for these parameters. In practice, we may need to choose a best one from the models mentioned in this paper to 7

balance accuracy and simplicity of calculations, and some relationships and properties of these models are also presented. The rest of the paper is organized as follows. In Section 2, basic assumptions are introduced and the main time redundant model is developed with three models as special cases. In Section 3, analysis of the main model and the three special cases are given, including number distribution of mission repeated tries, mission success probability, mean time to mission success and the optimal analysis for cost and revenue. In Section 4, properties and relationships of the main model and its three special cases are discussed for better understanding of them. In Section 5, some extended discussions are investigated briefly. In Section 6, numerical examples are provided to illustrate results given in the previous sections. Finally, conclusions and related discussions are summarized in Section 7.

2.

Model assumptions Assume that the time needed to finish a mission cannot be negligible, and the mission

window is a given interval [0,τ ] , which is long enough for the mission to be tried consecutively several times. That is to say, within the given mission window [0,τ ] , if the mission fails at the first time, then the mission can be tried again, if the mission fails again at the second time, then the mission can be tried at the third time, and so forth until it succeeds or stops at the end of the mission window. The mission may succeed or fail by time τ . Let

X i (i = 1, 2,…) be mission repeated durations (including mission failure durations or the last mission success duration or both), and they are independent and identically distributed random variables whose probability density function and distribution function are f (t ) and

F (t ) , respectively. Assume that X i for i = 1, 2,… have a common lower bound 0 and a common upper bound

b , which means

f (t ) = 0, F (t ) = 0

when

t<0

and

f (t ) = 0, F (t ) = 1 when t > b . Then if a new mission repeated try can be conducted immediately after a failed one, there will be at least [τ / b] mission repeated tries in mission window [0,τ ] before the mission is failed by time τ ( [τ / b] means the largest integer 8

less than or equal to τ / b ). In practice, sometimes τ can be determined by letting [τ / b] reach some threshold so that enough mission repeated tries can be carried out in this mission window. The main model of this paper is proposed in the following context together with three of its special cases. Before that, for simplicity, we denote the number of mission repeated tries in interval [0, t ] as N (t ) for the main model. It is obvious that N (t ) is a counting process determined by distribution function F (t ) and other assumptions. Note that N (t ) is a constant for t ≥ τ because the mission repeated tries will stop at the mission success time or at the end of the given mission window. Besides, N (t ) is denoted by N j (t ) in Model j

( j = 1, 2,3) . Main Model. The ith mission repeated try fails definitely if its duration is no greater than a given value a (short mission repeated try) and fails with probability qki = r i q, 0 < r ≤ 1 k

if its duration is greater than a (long mission repeated try), where ki is the number of mission repeated tries before the ith one whose durations are greater than a . The ith repair duration is denoted as Yi and assume that Yi

(i = 1, 2,…) are independent and

identically distributed random variables with density function g (t ) and distribution function

G (t ) . With consideration of the counting process N (t ) , the case that the mission has succeeded by time τ is shown in Figure 2 and the case that the mission has failed by time

τ is shown in Figure 3. (Figure 2 insert here about) (Figure 3 insert here about) In Figure 2, the first mission repeated try has a duration greater than a and it fails with probability q ; The second one has a duration less than a and it fails definitely; The third one has a duration greater than a and it fails with probability rq ; Finally, the nth one has a duration greater than a and it succeeds with probability 1 − r n q before time τ . k

Figure 3 is similar to Figure 2, except that the nth 9

mission repeated try fails with

probability r n q and the (n + 1)th one does not end before time τ . k

Remark 1. Three special cases of the main model are given as follows: Model 1. In the main model, let r = 1 , which means the failure probability qki of each long mission repeated try is a constant q . Then the case that the mission is successful by time τ is shown in Figure 4. Note that the case that the mission is failed by time τ can be given similarly thus omitted for the three models. Figure 4 is similar to Figure 2, except that the

nth mission repeated try succeeds with probability p rather than 1 − r kn q . (Figure 4 insert here about) Model 2. In Model 1, let Yi = 0 , which implies that repair of each failure can be done instantly without consuming time resources. Then for Model 2, the case that the mission is successful by time τ is shown in Figure 5. Figure 5 is similar to Figure 4, except that there is no mission repair duration in it. (Figure 5 insert here about) Model 3. In Model 2, let a = 0 , which indicates that all mission repeated tries are assumed to fail with probability q rather than 1 . Then the case that the mission is successful by time

τ is shown in Figure 6 for Model 3. Figure 6 is similar to Figure 5, except that the second mission repeated try fails with probability q rather than 1 . (Figure 6 insert here about) Remark 2. Assumptions of these models are summarized as follows in Table 1. (Table 1 insert here about)

3. Analysis of the main model and its three special cases 3.1. Number of mission repeated tries Denote ϕ (t , n) = P{N (t ) = n} (t ≤ τ )

as the probability that there are exactly n

mission repeated tries in interval [0, t ] for the main model. Then ϕ (t , n) can be given as follows for the main model and its three special cases respectively. 3.1.1. Main model 10

For the main model, if there are n (n ≥ 1) mission repeated tries in interval [0, t ] , the first n − 1 mission repeated tries should all fail. For the case that the nth mission repeated try succeeds in interval [0, t ] (see Figure 1), we have,

PS {N (t ) = n} = pP{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t , X n > a, X 1 ≤ a, X 2 ≤ a, … , X n −1 ≤ a} n

k −1

k =2

l =1

+ ∑ (∏ ql −1 ) pk −1 P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t , X n > a,k − 1 of X 1 , X 2 , … , X n −1 are greater than a} = pP{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t , X n > a, X 1 ≤ a, X 2 ≤ a, … , X n −1 ≤ a} k −1 n  n − 1 + ∑ [∏ (r l −1q )](1 − r k −1q )   P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t , X 1 , k = 2 l =1  k − 1 X 2 … , X n − k ≤ a, X n − k +1 , X n − k + 2 … , X n > a}

∫

=

P{a < X n ≤ t − (u1 + u2 + L + un −1 ) − (v1 + v2 +L + vn −1 )}

0
dP{ X 1 ≤ u1 , X 1 ≤ a}dP{ X 2 ≤ u2 , X 2 ≤ a}L dP{ X n −1 ≤ un −1 , X n −1 ≤ a} d {Y1 ≤ v1}L d {Yn −1 ≤ vn −1} n  n − 1 + ∑ r ( k −1)( k − 2)/2 q k −1 (1 − r k −1q )  P{a < X n ≤ t − (u1 + u2 + L  ∫ k =2  k − 1 0
dP{ X n − k ≤ un − k , X n − k ≤ a}dP{a < X n − k +1 ≤ un − k +1}L dP{a < X n −1 ≤ un −1}d {Y1 ≤ v1}L d {Yn −1 ≤ vn −1} n  n − 1 = pF n −1 (a ) F (a ) F1( n −1) ∗ F2 ∗ G ( n −1) (t ) + ∑ r ( k −1)( k − 2)/2 q k −1 (1 − r k −1q )   k =2  k − 1 × F n − k (a ) F k (a ) F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ) n  n − 1 n−k k ( n−k ) = ∑ r ( k −1)( k − 2)/2 q k −1 (1 − r k −1q )  ∗ F2( k ) ∗ G ( n −1) (t ).  F (a ) F (a ) F1 − k 1 k =1  

(3.1) where

p = 1− q

t

and

F ( n ) (t ) = ∫ F ( n −1) (t − u )dF (u ) 0

distribution function F (t ) with itself. Note that F

(0)

is a

n -fold convolution of

(t ) = 1 and F (1) (t ) = F (t ) for t ≥ 0 .

In (3.1), the first equality comes from the fact that there can be k (k = 1, 2, … , n) mission repeated tries whose durations are greater than a when N (t ) = n and the last mission succeeds, and the only long mission repeated try must be the last one for the case that k = 1 . For the case that the nth mission repeated try fails, we have, 11

PF {N (t ) = n} = P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t < X 1 + Y1 + L X n + Yn + X n +1 , X 1 ≤ a, n

k

k =1

l =1

X 2 ≤ a, … , X n ≤ a} + ∑ (∏ ql −1 ) P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t < X 1 + Y1 + L X n + Yn + X n +1 , k of X 1 , X 2 , … , X n are greater than a} = P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t < X 1 + Y1 + L X n + Yn + X n +1 , X 1 ≤ a, X 2 ≤ a,… , X n ≤ a} k n n + ∑ [∏ (r l −1q)]   P{ X 1 + Y1 + L X n −1 + Yn −1 + X n ≤ t < X 1 + Y1 + L X n + Yn + X n +1 , k =1 l =1 k  X 1 , X 2 … , X n − k ≤ a, X n − k +1 , X n − k + 2 … , X n > a}

∫

=

P{Yn + X n +1 > t − (u1 + u2 + L + un ) − (v1 + v2 + L + vn −1 )}

0
dP{ X 1 ≤ u1 , X 1 ≤ a}L dP{ X n ≤ un , X n ≤ a}d {Y1 ≤ v1}L d {Yn −1 ≤ vn −1} n n +∑ r k ( k −1)/2 q k   P{Yn + X n +1 > t − (u1 + u2 + L + un ) − (v1 + v2 + L ∫ k =1  k  0
+ vn −1 )}dP{ X 1 ≤ u1 , X 1 ≤ a}L dP{ X n − k ≤ un − k , X n − k ≤ a}dP{a <

X n − k +1 ≤ un − k +1}L dP{a < X n ≤ un }d{Y1 ≤ v1}L d{Yn −1 ≤ vn −1} n n = F n (a) ∫ [1 − G ∗ F (t − u )]dF1( n ) ∗ G ( n −1) (u )+∑ r k ( k −1)/2 q k   F n − k (a) F k (a) k =1 k  0 t

t

× ∫ [1 − G ∗ F (t − u )]dF1( n − k ) ∗ F2( k ) ∗ G ( n −1) (u ) 0 n

= ∑r

t  n  n−k k q   F (a) F (a) ∫ [1 − G ∗ F (t − u )]dF1( n − k ) ∗ F2( k ) ∗ G ( n −1) (u ). k  0

k ( k −1)/2 k

k =0

(3.2) In (3.2), the first equality comes from the fact that there can be k (k = 0,1, … , n) mission repeated tries whose durations are greater than a when N (t ) = n and the last mission fails, and all the mission repeated tries must be short ones for the case that k = 0 . Then the probability that there are n mission repeated tries in interval [0, t ] can be given as follows,

P{N (t ) = n} = PS {N (t ) = n} + PF {N (t ) = n} n  n − 1 n−k k ( n−k ) = ∑ r ( k −1)( k − 2)/2 q k −1 (1 − r k −1q )  ∗ F2( k ) ∗ G ( n −1) (t ) (3.3)  F (a ) F (a ) F1 k =1  k − 1 t n n + ∑ r k ( k −1)/2 q k   F n − k (a ) F k (a ) ∫ [1 − G ∗ F (t − u )]dF1( n − k ) ∗ F2( k ) ∗ G ( n −1) (u ). k =0 k  0

Then ϕ (t , n) = P{N (t ) = n} can be given by its Laplace-Stieltjes transform as follows, 12

 n − 1 ( k −1)( k − 2)/2 (1 − r k −1q )[ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k r k − 1 k =1   n

ϕˆ ( s, n) = Gˆ n −1 ( s ) F (a) Fˆ2 ( s )∑ 

n n + Gˆ n −1 ( s )[1 − Fˆ ( s )Gˆ ( s )]∑   r k ( k −1)/2 [ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k . k =0  k 

(3.4) Note that ϕˆ ( s, 0) = P{N (t ) = 0} = P{ X 1 > t} = F (t ) .

(3.5)

2

Expectations of N (t ) and N (t ) can be given by their Laplace-Stieltjes transforms as follows,

Eˆ [ N ( s )] ∞ n  n − 1 ( k −1)( k − 2)/2 = F (a ) Fˆ2 ( s )∑ nGˆ n −1 ( s )∑  (1 − r k −1q )[ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k −1 r n =1 k =1  k − 1  ∞ n n + [1 − Fˆ ( s )Gˆ ( s )]∑ nGˆ n −1 ( s )∑   r k ( k −1)/2 [ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k , n =1 k =0  k 

(3.6)

Eˆ [ N 2 ( s )] ∞ n  n − 1 ( k −1)( k − 2)/2 (1 − r k −1q )[ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k −1 = F (a ) Fˆ2 ( s )∑ n 2Gˆ n −1 ( s )∑  r n =1 k =1  k − 1  ∞ n n + [1 − Fˆ ( s )Gˆ ( s )]∑ n 2Gˆ n −1 ( s )∑   r k ( k −1)/2 [ F (a ) Fˆ1 ( s )]n − k [qF (a ) Fˆ2 ( s )]k . n =1 k =0  k 

(3.7) Remark 3. Variance of N (t ) can be given easily by its Laplace-Stieltjes transform as

Vˆ [ N ( s )] = Eˆ [ N 2 ( s )] − Eˆ 2 [ N ( s )] . Higher-order moments of N (t ) can be given in a similar way as (3.6) and (3.7), thus omitted here.

3.1.2. Model 1 Let

r =1

in

the

main

model,

then

for

Model

1,

the

probability

ϕ1 (t , n) = P{N1 (t ) = n} that there are n mission repeated tries in interval [0, t ] can be given as follows,

P{N1 (t ) = n} = PS {N1 (t ) = n} + PF {N1 (t ) = n}, 13

(3.8)

where n  n − 1 n−k k ( n−k ) PS {N1 (t ) = n} = ∑ q k −1 p  ∗ F2( k ) ∗ G ( n −1) (t ),  F (a ) F (a ) F1 k − 1 k =1   t n n PF {N1 (t ) = n} = ∑ q k   F n − k (a ) F k (a ) ∫ [1 − G ∗ F (t − u )]dF1( n − k ) ∗ F2( k ) ∗ G ( n −1) (u ), k =0 k  0

(3.9) correspond to the case that the nth mission repeated try succeeds and fails, respectively. Then Laplace-Stieltjes transform of ϕ1 (t , n) = P{N1 (t ) = n} can be presented as follows,

 n − 1 k −1 n − k k ˆ n −1 ˆ n−k ˆk  q pF (a ) F (a ) F1 ( s ) F2 ( s )G ( s ) k =1  k − 1  n

ϕˆ1 ( s, n) = ∑ 

n n + ∑   q k F n − k (a ) F k (a ) Fˆ1n − k ( s ) Fˆ2k ( s )Gˆ n −1 ( s ) k =0  k  n n − ∑   q k F n − k (a ) F k (a ) Fˆ1n − k ( s ) Fˆ2k ( s ) Fˆ ( s )Gˆ n ( s ) k =0  k 

(3.10)

= pGˆ n −1 ( s ) F (a ) Fˆ2 ( s )[ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]n −1 + [1 − Fˆ ( s )Gˆ ( s )]Gˆ n −1 ( s )[ F (a ) Fˆ ( s ) + qF (a ) Fˆ ( s )]n . 1

2

That means,

ϕˆ1 ( s, n) = Gˆ ( s )[ F (a) Fˆ1 ( s ) + qF (a) Fˆ2 ( s )]ϕˆ1 ( s, n − 1),

n ≥ 2, ϕˆ1 ( s,1) = pF (a) Fˆ2 ( s ) + [1 − Fˆ ( s )Gˆ ( s )][ F (a) Fˆ1 ( s ) + qF (a) Fˆ2 ( s )]. Note that ϕˆ ( s, 0) = P{N (t ) = 0} = P{ X 1 > t} = F (t ) .

(3.11)

(3.12)

2

Expectations of N1 (t ) and N1 (t ) can be given by their Laplace-Stieltjes transforms as follows, ∞

Eˆ [ N1 ( s )] = ∑ nϕˆ1 ( s, n) n =1 ∞

= ∑ n{Gˆ ( s )[ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]}n −1ϕˆ1 ( s,1) n =1

= { pF (a ) Fˆ2 ( s ) + [1 − Fˆ ( s )Gˆ ( s )][ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]} × [1 − F (a ) Fˆ ( s )Gˆ ( s ) − qF (a ) Fˆ ( s )Gˆ ( s )]−2 , 1

2

and

14

(3.13)

∞

Eˆ [ N12 ( s )] = ∑ n 2ϕˆ1 ( s, n) n =1 ∞

= ∑ n 2 [1 + F (a ) Fˆ1 ( s )Gˆ ( s ) + qF (a ) Fˆ2 ( s )Gˆ ( s )]n −1ϕˆ1 ( s,1) n =1

= { pF (a ) Fˆ2 ( s ) + [1 − Fˆ ( s )Gˆ ( s )][ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]} ×[1 + F (a ) Fˆ ( s )Gˆ ( s ) + qF (a ) Fˆ ( s )Gˆ ( s )] 1

(3.14)

2

×[1 − F (a ) Fˆ1 ( s )Gˆ ( s ) − qF (a ) Fˆ2 ( s )Gˆ ( s )]−3 . 3.3.3. Model 2 Let Yi = 0 in Model 1, which means Gˆ ( s ) = 1 , then for Model 2, the probability

ϕ2 (t , n) = P{N 2 (t ) = n} that there are n mission repeated tries in interval [0, t ] can be given as follows,

P{N 2 (t ) = n} = PS {N 2 (t ) = n} + PF {N 2 (t ) = n},

(3.15)

where n  n − 1 n−k ( n−k ) k PS {N 2 (t ) = n} = ∑ q k −1 p  ∗ F2( k ) ,  F (a ) F (a ) F1 k − 1 k =1   (3.16) t n ( n−k ) (k ) k n n−k k PF {N 2 (t ) = n} = ∑ q   F (a ) F (a ) ∫ [1 − F (t − u )]dF1 ∗ F2 (u ), k =0 k  0

correspond to the case that the nth mission repeated try succeeds and fails, respectively. Then Laplace-Stieltjes transform of ϕ 2 (t , n) = P{N 2 (t ) = n} can be presented as follows,

ϕˆ2 ( s, n)= pF (a) Fˆ2 ( s )[ F (a) Fˆ1 ( s ) + qF (a) Fˆ2 ( s )]n −1 + [1 − Fˆ ( s )][ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]n .

(3.17)

That means,

ϕˆ2 ( s, n) = [ F (a) Fˆ1 ( s ) + qF (a) Fˆ2 ( s )]ϕˆ2 ( s, n − 1),

n ≥ 2, ϕˆ2 ( s,1) = pF (a) Fˆ2 ( s ) + [1 − Fˆ ( s )][ F (a) Fˆ1 ( s ) + qF (a) Fˆ2 ( s )]. Note that ϕˆ2 ( s, 0) = P{N 2 (t ) = 0} = P{ X 1 > t} = F (t ) .

(3.18)

(3.19)

2

Expectations of N 2 (t ) and N 2 (t ) can be given by their Laplace-Stieltjes transforms as follows,

Eˆ [ N 2 ( s )] = { pF (a ) Fˆ2 ( s ) + [1 − Fˆ ( s )][ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]} × [1 − F (a ) Fˆ ( s ) − qF (a ) Fˆ ( s )]−2 , 1

2

15

(3.20)

Eˆ [ N 22 ( s )] = { pF (a ) Fˆ2 ( s ) + [1 − Fˆ ( s )][ F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )]} ×[1 + F (a ) Fˆ1 ( s ) + qF (a ) Fˆ2 ( s )][1 − F (a ) Fˆ1 ( s ) − qF (a ) Fˆ2 ( s )]−3 .

(3.21)

3.1.4 Model 3 Let a = 0 in Model 2, then for Model 3, the probability ϕ (t , n) = P{N 3 (t ) = n} that there are n mission repeated tries in interval [0, t ] can be given as follows,

P{N 3 (t ) = n} = PS {N 3 (t ) = n} + PF {N 3 (t ) = n} t

= q n −1 pF ( n ) (t ) + ∫ q n [1 − F (t − u )]dF ( n ) (u ) 0

=q

n −1

n −1

pF

=q F

(n)

(n)

(t ) + q [ F n

(t ) − q F n

(n)

( n +1)

(t ) − F

( n +1)

(3.22)

(t )]

(t ),

where

PS {N 3 (t ) = n} = q n −1 pP{ X 1 + X 2 + L + X n ≤ t} = q n −1 pF ( n ) (t ), PF {N 3 (t ) = n} = q n P{ X 1 + X 2 + L + X n ≤ t < X 1 + X 2 + L + X n +1} t

= q n ∫ P{ X n +1 > t − u}dP{ X 1 + X 2 + L + X n ≤ u}

(3.23)

0 t

= q ∫ [1 − F (t − u )]dF ( n ) (u ), n

0

correspond to the case that the nth mission repeated try succeeds and fails, respectively. Then Laplace-Stieltjes transform of ϕ3 (t , n) = P{N 3 (t ) = n} can be presented as follows,

ϕˆ3 ( s, n)= q n −1 Fˆ n ( s) − q n Fˆ n +1 ( s) = [1 − qFˆ ( s)]q n −1 Fˆ n ( s). Note that ϕˆ3 ( s, 0) = P{N 3 (t ) = 0} = P{ X 1 > t} = F (t ) .

(3.24) (3.25)

2

Expectations of N 3 (t ) and N 3 (t ) can be given by their Laplace-Stieltjes transforms as follows, ∞

Eˆ [ N 3 ( s )] = ∑ nϕˆ3 ( s, n) n =1

∞

= [1 − qFˆ ( s )]Fˆ ( s )∑ n[qFˆ ( s )]n −1 n =1

= Fˆ ( s )[1 − qFˆ ( s )]−1 ,

16

(3.26)

∞

Eˆ [ N 32 ( s )] = ∑ n 2ϕˆ3 ( s, n) n =1

∞

= [1 − qFˆ ( s )]Fˆ ( s )∑ n 2 [qFˆ ( s )]n −1

(3.27)

n =1

= Fˆ ( s )[1 + qFˆ ( s )][1 − qFˆ ( s )]−2 .

3.2. Mission success probability Denote the probability that the mission succeeds before time t (t ≤ τ ) as H (t ) for the main model, and we have H (t ) = H j (t ) for Model j ( j = 1, 2,3) . Since there can be

n (n = 1, 2,… , ∞) mission repeated tries before time t , from Section 3.1, we have, ∞

∞

n =1

n =1

H (t ) = ∑ PS {N (t ) = n} = ∑ Pr{N (t ) = n and the nth mission repeated try successes}, (3.28) which means, ∞ n  n − 1 ( k −1)( k − 2)/2 k −1 H (t ) = ∑∑  q (1 − r k −1q ) F n − k (a ) F k (a ) F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ), r n =1 k =1  k − 1 

(3.29) and ∞ n  n − 1 k −1 n − k ( n−k ) k H1 (t ) = ∑∑  ∗ F2( k ) ∗ G ( n −1) (t ),  q pF (a ) F (a ) F1 n =1 k =1  k − 1  ∞ n  n − 1 k −1 n − k ( n−k ) k H 2 (t ) = ∑∑  ∗ F2( k ) (t ),  q pF (a ) F (a ) F1 n =1 k =1  k − 1 

(3.30)

∞

H 3 (t ) = ∑ q n −1 pF ( n ) (t ). n =1

Denote the time to mission success as T for the main model and

T = T j for Model

j ( j = 1, 2,3) , then their k -order moment E (T k ) , k = 1, 2,… can be given as follows, ∞

dk E (T ) = ∫ t dH (t ) = (−1) ds k 0 k

k

k

where for the main model, 17

s =0

Hˆ ( s ),

(3.31)

∞ n  n − 1 ( k −1)( k − 2)/2 k −1 Hˆ ( s ) = ∑∑  q (1 − r k −1q )[ F (a ) Fˆ1 ( s )]n − k [ F (a ) Fˆ2 ( s )]k Gˆ n −1 ( s ) r k − 1 n =1 k =1   (3.32) ∞ n  n − 1 ( k −1)( k − 2)/2 k −1 k −1 n − k k n − 1 = ∑∑  q (1 − r q )[ F (a ) Fˆ1 ( s )] [ F (a ) Fˆ2 ( s )] Gˆ ( s ). r n =1 k =1  k − 1 

And for Model j ( j = 1, 2,3) , ∞ n  n − 1 k −1 n−k k ˆ n −1 ˆ ˆ Hˆ 1 ( s ) = ∑∑   q p[ F (a ) F1 ( s )] [ F (a ) F2 ( s )] G ( s ) − 1 k n =1 k =1   = pF (a ) Fˆ ( s )[1 − F (a ) Fˆ ( s )Gˆ ( s ) − qF (a ) Fˆ ( s )Gˆ ( s )]−1 , 2

1

2

 n − 1 k −1 n−k k ˆ ˆ Hˆ 2 ( s ) = ∑∑   q p[ F (a ) F1 ( s )] [ F (a ) F2 ( s )] − k 1 n =1 k =1   = pF (a ) Fˆ ( s )[1 − F (a ) Fˆ ( s ) − qF (a ) Fˆ ( s )]−1 , ∞

n

1

2

(3.33)

2

∞

Hˆ 3 ( s ) = ∑ q n −1 pFˆ ( n ) ( s ) = pFˆ ( s )[1 − qFˆ ( s )]−1. n =1

In practice, parameters in these models are determined in the process of system design with the help of data analysis, which means we can provide some guidance in that process by discussing the impact of these parameters on the mission success probability. Obviously,

H (t ) is a function of time t and in the following context, we will illustrate that it is an increasing function of t and a decreasing function of parameters q, r . That means the mission success probability for the main model will increase if we are given a longer mission window or if we can lower the failure probability of each mission repeated try. Note that for distribution functions F1 (t ) and F2 (t ) , their convolution can be given as t

F1 ∗ F2 (t ) = ∫ F1 (t − u )dF2 (u )

and for density functions

f1 (t )

and

f 2 (t ) , their

0 t

convolution can be given as f1 ∗ f 2 (t ) = ∫ f1 (t − u ) f 2 (u )du . 0

For the main model, H (t ) is an increasing function of t because its derivative function with respect to t is nonnegative for all t ≥ 0 as follows, ∞ n  n − 1 ( k −1)( k − 2)/2 k −1 h(t ) = ∑∑  q (1 − r k −1q ) F n − k (a ) F k (a ) f1( n − k ) ∗ f 2( k ) ∗ g ( n −1) (t ). (3.34) r − k 1 n =1 k =1  

We can also find that H (t ) is a decreasing function of parameter q because its partial 18

derivative function L(t ) =

∂H (t ) with respect to q is nonpositive for all t ≥ 0 . Related ∂q

discussions can be given as follows: ∞ n  n − 1 ( k −1)( k − 2)/2 k −1 H (t ) = ∑∑  q (1 − r k −1q ) F n − k (a ) F k (a ) F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ) r n =1 k =1  k − 1  ∞ ∞  n − 1 ( k −1)( k − 2)/2 k −1 = ∑∑  q (1 − r k −1q ) F n − k (a ) F k (a ) F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ), r − k 1 k =1 n = k  

(3.35) which means, ∞ ∞  n − 1 n−k k L(t ) = ∑ [r ( k −1)( k − 2)/2 (k − 1)q k − 2 − r ( k −1) k /2 kq k −1 ]∑   F (a) F (a) k =2 n = k  k − 1 × F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ).

(3.36)

Taking the Laplace-Stieltjes transform of L(t ) , we have, ∞ ∞  n − 1 n−k ˆ Lˆ ( s ) = ∑ [r ( k −1)( k − 2)/2 (k − 1)q k − 2 − r ( k −1) k /2 kq k −1 ]∑   [ F (a ) F1 ( s )] k =1 n = k  k − 1 ×[ F (a ) Fˆ ( s )]k G n −1 ( s ) 2

∞

= ∑ [r ( k −1)( k − 2)/2 (k − 1)q k − 2 − r ( k −1) k /2 kq k −1 ] k =1

[ F (a ) Fˆ2 ( s )]k G k −1 ( s ) . [1 − F (a ) Fˆ1 ( s )G ( s )]k (3.37)

Then by taking the inverse Laplace-Stieltjes transform, we have, ∞

L(t ) = ∑ [r ( k −1)( k − 2)/2 (k − 1)q k − 2 − r ( k −1) k /2 kq k −1 ]Q2( k ) ∗ G ( k −1) (t ) − Q2 (t ) k =2

∞

∞

k =2

k =2

≤ ∑ r ( k −1)( k − 2)/2 (k − 1)q k − 2Q2( k −1) ∗ G ( k − 2) (t ) − ∑ r ( k −1) k /2 kq k −1Q2( k ) ∗ G ( k −1) (t ) − Q2 (t ) = Q2 (t ) − Q2 (t ) = 0. (3.38) ∞

where

Q2 (t ) = ∑ F n −1 (a ) F (a ) F1( n −1) ∗ F2 ∗ G ( n −1) (t )

has

Laplace

transform

n =1

∞

Qˆ 2 ( s ) = F (a ) F2 ( s )∑ F n −1 (a ) Fˆ1( n −1) ( s )Gˆ n −1 ( s ) = n =1

F (a ) Fˆ2 ( s ) . Obviously Q2 (t ) is 1 − F (a ) Fˆ1 ( s )G ( s )

a distribution function, because Q2 (t ) = 0 , Q2 (∞) =

∞

∑F n =1

19

n −1

(a) F (a) =

F (a) = 1 and 1 − F (a)

∞

Q2 (t ) has a nonnegative derivative function q2 (t ) = ∑ F n −1 (a ) F (a ) f1( n −1) ∗ f 2 ∗ g ( n −1) (t ) . n =1

Note that the inequality is true because P{ X + Y ≤ t} must be less than or equal to

P{ X ≤ t}

for

any

positive

random

variable

Y,

which

implies

Q2( k ) ∗ G ( k −1) (t ) ≤ Q2( k −1) ∗ G ( k − 2) (t ) . Similarly, we can also find that H (t ) is a decreasing function of parameter r because its partial derivative function J (t ) =

∂H (t ) with respect to r is non-positive for all t ≥ 0 . ∂r

Similar to (3.35)-(3.36), we get ∞ ∞  n − 1 J (t ) = ∑ [(k − 1)(k − 2)r ( k −1)( k − 2)/2−1q k −1 − (k − 1)kr ( k −1) k /2−1q k ]∑   k =2 n = k  k − 1 × F n − k (a ) F k (a ) F1( n − k ) ∗ F2( k ) ∗ G ( n −1) (t ).

(3.39)

Taking the Laplace-Stieltjes transform of J (t ) , we have, ∞ ∞  n − 1 Jˆ ( s ) = ∑ [(k − 1)(k − 2)r ( k −1)( k − 2)/2−1q k −1 − (k − 1)kr ( k −1) k /2−1q k ]∑   k =1 n = k  k − 1 × [ F (a ) Fˆ ( s )]n − k [ F (a ) Fˆ ( s )]k G n −1 ( s ) 1

2

∞

= ∑ [(k − 1)(k − 2)r ( k −1)( k − 2)/2−1q k −1 − (k − 1)kr ( k −1) k /2−1q k ] k =1

[ F (a ) Fˆ2 ( s )]k G k −1 ( s ) . [1 − F (a ) Fˆ1 ( s )G ( s )]k (3.40)

Then by taking the inverse Laplace-Stieltjes transform, we have, ∞

J (t ) = ∑ [(k − 1)(k − 2)r ( k −1)( k − 2)/2−1q k −1 − (k − 1)kr ( k −1) k /2−1q k ]Q2( k ) ∗ G ( k −1) (t ) k =1

∞

∞

k =1

k =2

≤ ∑ (k − 1)(k − 2)r ( k −1)( k − 2)/2−1q k −1Q2( k −1) ∗ G ( k − 2) (t ) − ∑ (k − 1)kr ( k −1) k /2−1q k Q2( k ) ∗ G ( k −1) (t ) = 0. (3.41) where Q2 (t ) is the same to (3.38). Obviously, these properties can be applied to Models 1, 2 and 3 because they are all special cases of the main model. 20

4. Relationships of the main model and its three special cases In this section, we will simply discuss properties of the main model and the three special cases, and make some comparisons which may help us understand them better. Relationships and properties of the models are summarized as follows. (1) When r = 1 , the main model will reduce to Model 1. When Yi = 0 for all i = 1, 2,… , Model 1 will reduce to Model 2. When a = 0 , Model 2 will reduce to Model 3. That is to say,

Model 3 ⊆ Model 2 ⊆ Model 1 ⊆ main model . (2) N (0) = 0 , and N (t ) will be constants when t ≥ τ because the mission repeated tries will stop at time τ or before time τ . (3) T ≤ T1 and N1 (t ) ≥ N (t ) with probability 1 , because the main model has greater mission success probabilities pki = 1 − r i q ≥ p (0 < r < 1) than Model 1 which makes k

it stop earlier than Model 1 with probability 1 . (4) T2 ≤ T1 and

N 2 (t ) ≥ N1 (t ) with probability 1 , because Model 1 needs an extra

duration Yi ≥ 0 to start another mission repeated try after the ith mission failed duration, which leads to a later mission success time than Model 2. (5) T3 ≤ T2 and N 2 (t ) ≥ N 3 (t ) with probability 1 , because obviously mission success time of Model 3 must be earlier than the mission success time of Model 2, and at any time

N 3 (t ) increases 1 , N 2 (t ) also increases 1 . (6) From (3)-(5), we have, E[ N 2 (t )] ≥ E[ N1 (t )] ≥ E[ N (t )] , E[ N 2 (t )] ≥ E[ N 3 (t )] and

E (T ) ≤ E (T1 ) , E (T3 ) ≤ E (T2 ) ≤ E (T1 ) . (7) H (0) = 0 , H (∞) = 1 , and H (t ) increases with t ≥ 0 . For a > 0 , H (t ) = 0 when

0
21

means H1 (t ) = P{T1 ≤ t} ≤ P{T ≤ t} = H (t ) . Then similarly from (4)-(5), we also have,

H1 (t ) ≤ H 2 (t ) ≤ H 3 (t ) for t ≥ 0 . The properties above show that Models 1-3 are given by simplifying the main model to different extents. In practice, we may need to choose one from four models to balance accuracy and simplicity of calculations, and those properties can provide some guidance we need for model selection.

5.

Some extended discussions of the main model

5.1. The Laplace transform method for the main model In Section 3, Laplace transform and inverse Laplace transform should be used to get measures such as the number of mission repeated tries and the mission success probability. For numerical calculation of the main model and its three special cases, Maple software can be used to get those Laplace transforms and inverse Laplace transforms. As mentioned in Section 1, the distributions of mission repeated durations and repair durations are given by fitting real data in practice, therefore we should choose those distributions that have analytical forms of Laplace transforms if possible, and get these Laplace transforms by using simple built-in functions in Maple. However, after several steps of calculation, the results in the last step are always complex and do not have analytic inverse Laplace transforms. Then some numerical methods for Laplace transform inversion should be employed (Cohen, 2007). Note that in Section 3, formulas are given by their Laplace-Stieltjes transforms rather than Laplace transforms. The Laplace-Stieltjes transform of a function F (t ), t ≥ 0 is defined as ∞

∞

Fˆ ( s ) = ∫ e − st dF (t ) and its relationship with the Laplace transform F ∗ ( s ) = ∫ e − st F (t )dt 0

0

* can be given by Fˆ ( s ) = sF ( s ) − F (0) (Kuhfittig, 2013).

In this paper, we use the numerical Laplace transform inversion formula proposed by Stehfest (1970b) rather than a wrong version in Stehfest (1970a), which is a very efficient formula given as follows,

22

ln 2 2 M  ln 2 F (t ) = V j ,M F *  ∑ t j =1  t

 j , V j , M = (−1) M + j 

min( j , M )

∑

k =[

j +1 ] 2

k M (2k )! , ( M − k )!k !(k − 1)!( j − k )!(2k − j )! (5.1)

is often empirically chosen from 4,5, 6, 7,8 and here we use M = 7 after

where M

some trials. This Laplace inversion method, which is called the Gaver-Stehfest algotithm, has been widely used for its simplicity and good performance, and its convergence problem has been concluded and discussed in Kuznetsov (2013).

5.2. The Monte Carlo method for the extended model In practice, there may be different kinds of failure modes for a mission, which indicates that the main model can be extended to a more general model (called the extended model). In this section, we assume that the ith mission repeated duration X% i (i = 1, 2,… , ∞) can be determined by the first one of m different kinds of failures. That implies the failure mode will be ei where function

Z vi

(ei = 1, 2,… , m) and X% i will be Z ei ,i if Z ei ,i = min( Z1i , Z 2i ,… , Z mi ) ,

(v = 1, 2,… , m) has probability density function

F%v (t ) . The ith

distribution function

repair duration

G% ei (t ) if

X% i = Z ei ,i

Y%i

f%v (t ) and distribution

has probability function

g% ei (t ) and

(ei = 1, 2,… , m) . Note that

0 ≤ Z vi ≤ b

(v = 1, 2,… , m) because it is assumed that 0 ≤ X% i ≤ b . We also assume that the same kind of long failures have decreasing failure probability, that means, the ith mission repeated try fails definitely if its duration is no greater than a given value a and otherwise fails with k% probability q%ei ( i ) = r%ei ei ( i ) q%ei (0 < r%ei ≤ 1, ei = 1, 2,… , m) , where k%v (i ) (v = 1, 2,… , m) is the

number of long mission repeated tries before the ith one that have failure mode v and

r%v , q%v are parameters corresponding to failure mode v

(v = 1, 2,… , m) . Obviously,

k%v (1) = 0 and k%v (1) ≤ k%v (2) ≤ L ∞ . To analyze this extended model, the Laplace transform method in Section 3 will be very complex for large m , then the Monte Carlo simulation method will be used in this section 23

instead. Main steps to calculate mean number of mission repeated tries E[ N% (τ )] and mission success probability H% (τ ) are given as follows: Step 1. Set initial value j = 0 , and input a, b,τ , N , r%1 , r%2 ,… , r%m , q%1 , q%2 ,… , q%m . Step 2. If j ≤ N , let i = 0 , T = 0 , k%1 = k%2 = L = k%m = 0 . If j = N + 1 , turn to Step 7. Step 3. Generate z1i , z2i ,… , zmi from distributions F%1 (t ), F%2 (t ),… , F%m (t ) , respectively. k% Step 4. For zei ,i = min( z1i , z2i ,… , zmi ) , let k%ei = k%ei + 1 , q%ei = r%ei ei q%ei , T = T + zei ,i . If

T > τ , let x j = i , y j = 0 , j = j + 1 and turn to Step 2. If T ≤ τ , turn to Step 5. Step 5. If zei ,i ≤ a , generate wi from distribution G% ei (t ) and let T = T + wi , i = i + 1 , and then turn to Step 3. If zei ,i > a , turn to Step 6. Step 6. Generate ui from standard uniform distribution U [0,1] and compare ui with q%ei . If ui ≤ q%ei , generate wi from distribution G% ei (t ) and let T = T + wi , i = i + 1 , and turn to Step 3. If ui > q%ei , let x j = i + 1 , y j = 1 , and turn to Step 2. Step 7. Let x =

1 N

N

∑ xj , y = j =1

1 N

N

∑y

j

, and output x , y , namely Monte Carlo estimation

j =1

of E[ N% (τ )] and H% (τ ) . Note that the extended model will be the main model if we set m = 1 .

5.3. Comparison of the Laplace transform method and the Monte Carlo method For calculation of the main model, the Laplace method is much more efficient than the Monte Carlo method because it needs much less time for calculation. This advantage becomes more obvious for larger τ , because calculation times of the Monte Carlo method depends on

τ while calculation times of the Laplace method will never change with τ . Comparing with the Monte Carlo method, the results given by the Laplace transform method are accuate enough, which will be showed by numerical examples in Section 6.1. However, the Laplace method does not always works, for distributions without analytical 24

form of Laplace transforms, approximation of their Laplace transforms can be used, which may cause calculation errors to some extent. Therefore, to use the model in this paper, we prefer to choose distributions that have analytical form of Laplace transforms to fit practical data, which is reasonable in practice. Besides, when mission repeated durations and repair durations follow continuous distributions, H (a ) and E (a ) should be zero because the probability that the first mission repeated try has a length a is zero. For the Laplace transform method, this property (

H (a ) = E (a ) = 0 ) will lose to some extent in the process of using (5.1) to calculate inverse Laplace transform, but that does not have much effect on the accuracy of calculation for

t ≥ a which we are really interested in. In Section 3, we have concluded that E[ N (t )] and H (t ) are increasing functions of t , and this property will lose to some extent in the Monte Carlo method because accuracy of this method is no more than 1/ N and E[ N (t )] , H (t ) tend to converge for large t . The time needed for the Monte Carlo method is an increasing function of N , which means we may have to choose the Laplace method if the much time will be consumed for the accuracy we need.

5.4. Optimal analysis for cost and revenue In Section 3, we have discussed the impact of parameters on the mission success probability in our models to provide some guidance on system design. For the main model, the mission success probability H (τ ) is determined by the mission interval length τ , the mission repeated durations distribution F (t ) , the repair durations distribution G (t ) and parameters a, q, r . In practice, the mission interval and the shortest time needed for each mission repeated try can be assumed as fixed, which means τ and a are both constant. And the distributions F (t ) and G (t ) can also be determined according to historical data. Then to improve the mission success probability, measures should be conducted to minimize

q, r since H (τ ) is a decreasing function of parameters q and r . However, it is impossible to get lower failure probability in each mission repeated try without extra cost 25

invested for missions in practice, and the balancing of reliability measures and economical cost has been an important issue for a long time in the field of reliability (Chambari et al., 2012; Levitin et al., 2013; Qiu et al., 2017a; Qiu et al., 2017b). In this section, we further investigate optimization problems of the parameters with consideration of cost and revenue to provide a better guidance on system design in practice. We assume that the total cost in the fixed mission interval [0,τ ] includes two parts, the fixed maintenance cost for each mission repeated try and the extra cost for reducing the failure probability for each mission repeated try. Then the expected total cost in [0,τ ] can be represented as EC (τ ) = f1[ EN (τ )] + f 2 (q, r ) , where f1 is an increasing differentiable function of the expected number of mission repeated tries EN (τ ) and f 2 is an decreasing differentiable function of parameters q, r . Assume that the revenue in mission [0,τ ] is an increasing differentiable function of mission success probability H (τ ) , which is denoted as revenue in [0,τ ]

can be provided as

f 0 [ H (τ )] , then the expected net

ER(τ ) = f 0 [ H (τ )] − EC (τ ) . Now it becomes

important to decide the strategy of reducing the failure probability for each mission repeated try, namely to lower parameters q, r to a proper extent for a maximal net revenue ER (τ ) . The optimal problem can be presented as follows,

max ER(τ ) = f 0 [ H (τ )] − f1[ EN (τ )] − f 2 (q, r ) q ,r   s.t. 0 ≤ q ≤ 1, 0 ≤ r ≤ 1.

(5.2)

If the total funds for this mission are limited to c , a related constraint should be added to the optimal problem as follows,

max ER(τ ) = f 0 [ H (τ )] − f1[ EN (τ )] − f 2 (q, r )  q ,r  s.t. f1[ EN (τ )] − f 2 (q, r ) ≤ c,  0 ≤ q ≤ 1, 0 ≤ r ≤ 1. 

(5.3)

Note that with the increase of q or r , obviously the expected mission success revenue

f 0 [ H (τ )] , the expected mean maintenance cost f1[ EN (τ )] and the expected extra cost 26

f 2 (q, r ) will all decrease. When problem

can

always

be

solved

f 0 , f1 , f 2 are differentiable functions, this optimal since

point

( q, r )

lies

in

a

closed

set

S = {(q, r ) : f1[ EN (τ )] − f 2 (q, r ) ≤ c, 0 ≤ q ≤ 1, 0 ≤ r ≤ 1} . Obviously the optimal solution (q* , r * ) may be one point on boundaries of set S or a solution of the following series of equations in S ,

∂ ∂ ∂ ∂  ∂q ER(τ ) = ∂q f 0 [ H (τ )] − ∂q f1[ EN (τ )] − ∂q f 2 (q, r ),   ∂ ER(τ ) = ∂ f [ H (τ )] − ∂ f [ EN (τ )] − ∂ f (q, r ). 0 1 2  ∂r ∂r ∂r ∂r

(5.4)

To solve the optimal problem in this section, there are many heuristic search algorithms (like a genetic algorithm heuristic) that can be used (Levitin et al., 1994). For relatively simple

f 0 , f1 , f 2 , more information of this problem can be read from related plots. For example, as mentioned by Qiu et al. (2017b), the expected net revenue function some models can be defined as

R = a + b( A − Amin ) − E (C ) , where a is the fixed value of the contract,

b( A − Amin ) is the performance-based incentive and E (C ) . An application of the optimal analysis for cost and revenue can be found in the numerical examples in Section 6.

6.

Numerical examples In this section, we consider a TT&C system whose mission is transiting data from the

spaceflight to the observation station on the ground. Assume that the spaceflight passes overhead the observation in time interval [0,τ ] = [0,10] . The mission repeated durations

X i , i = 1, 2,…

follow

a

beta

distribution

with

probability

density

function

f ( x) = 0.8203125 x 2 (2 − x) 4 , 0 ≤ x ≤ 2 (here it means the parameter b = 2 ), and the repair durations Yi , i = 1, 2,… follow another beta distribution with probability density function g ( y ) = 20 y (1 − y ) , 0 ≤ y ≤ 1 . The shortest time needed for the mission to success 3

is assumed to be a = 1 , the failure probabilities for the long mission repeated tries whose 27

durations are greater than 1 are found to be q, rq, r q,… , where q = 0.7 and r = 0.5 2

for the main model. Then for the main model, we have,

E[ N (10)] ≈ 5.9373 and

H (10) = 0.6345 , which indicates that the TT&C is expected to carry out about 6 mission repeated tries in its mission interval [0,10] and the probability that the data transmission mission is successful at time 10 is 0.6345 .

6.1. Number of mission repeated tries and mission success probability To prove relationships of E[ N (t )] and H (t ) for the main model and its three special cases which are described by properties (6)-(8) in Section 4, their values for 0 ≤ t ≤ 100 are shown in Figure 7. Then it can be concluded that the number of mission repeated tries will greatly increase if repair of devices can be done instantly or failure probability for each long mission repeated try does not decrease. Besides, mission success probability of the TT&C system will increase if repair of devices can be done instantly or the failure probability for each long mission repeated try decreases quickly. Note that the system should stop at time

t = 10 , Figure 7 is given here for 0 ≤ t ≤ 100 only to see the convergence tends of

E[ N (t )] and H (t ) without being forced to stop. In Figure 7, When t is large enough, E[ N (t )] tends to be a constant and H (t ) tends to be one, which is reasonable because mission repeated durations X i (i = 1, 2,… , ∞) have a common upper bound b = 2 and the mission will succeed sooner or later if enough time is given. Take Model 3 as an example, the probability that the TT&C system is failed at time

t = 40 is less than

q t / b = 0.7 40/2 ≈ 0.008 , which means the mission success

probability at time t = 40 is no less than 0.992 . Note that N (t ) must stop increasing at the mission success time, which means E[ N (t )] also tends to convergent when t is large enough. (Figure 7 insert here about) To prove the correctness of the formulas in Section 3, results of E[ N (t )] and H (t ) 28

given by the Monte Carlo method (dotted lines) and the Laplace transform method (solid lines) are compared in Figure 8 for the main model and its three special cases. Since the Monte Carlo method is really time-assuming for large t , we only provide results for 0 ≤ t ≤ 20 and let N = 10000 . Note that the results by the Monte Carlo method may not be as accurate as the Laplace method since they have at most four digits when N = 10000 . But as shown in Figure 8, the values given by the two different methods are already very close, which verifies correctness of both methods. (Figure 8 insert here about)

6.2. Optimal analysis for cost and avenue Assume that the expected maintenance cost in interval [0,τ ] is f1[ EN (τ )] = 2 EN (τ ) , which means 2 million dollars are needed for each mission repeated try. The expected extra cost for reducing failure probabilities is f 2 (q, r ) = 10(1 − q ) + 0.5(r −1 − 1) , which is greater than

zero

for

q <1

or

r < 1 . The mean revenue in interval

[0,τ ]

is

f 0 [ H (τ )] = 5 H (τ ) + 15 , where 5 million dollars is for a successful mission and 15 million dollars is the fixed revenue obtained during the whole process. Then the net revenue in interval

[0,τ ]

will be

ER(τ ) = 5 H (τ ) − 2 EN (τ ) + 10q −0.5r −1 + 5.5 , where

0 ≤ q ≤ 1, 0.1 ≤ r ≤ 1 for practical reason. Note that costs for this mission are not enough and there is no limitation on the total expected cost. To consider the optimal parameter strategy for this system, values of ER (10) for different q, r are plotted in Figure 9 as follows. (Figure 9 insert here about) As shown in Figure 9. (left), ER (10) is an increasing function of q when r is small (for example,

r = 0.1, 0.2 ); and when r is large (for example, r = 0.3, 0.4,… ,1 ),

ER(10) increases first and then decreases with the increasing of q . Note that the decreasing of ER (10) for r = 0.3 cannot be observed easily but can be read from the data.

29

As shown in Figure 9. (right), ER (10) is an increasing function of r when q is small (for example, q = 0, 0.1, 0.2 ); and when q is large (for example, q = 0.3, 0.4,… ,1 ),

ER(10) increases first and then decreases with the increasing of q . Note that the decreasing of ER (10) for r = 0.3, 0.4 cannot be observed easily but can be read from the data. We can also conclude from Figure 9 that maximal ER (10) can be obtained when

r = 1 , which means failure probabilities are the same for all the mission repeated tries. Possible reason for this conclusion may be that the cost is relative high for upgrading device after each failure. Denote ER (10) as g (q ) for r = 1 and obviously g (q ) is an unimodal function of

q . Notice that

g (2) ≈ 1.7427, g (2.5) ≈ 1.7551, g (3) ≈ 1.7389 ,

which means the optimal solution of q lies in interval (0.2, 0.3) . By further calculation, we can shrink the interval to

(0.23, 0.25) , then to

(0.244, 0.246) , and to

(0.2448, 0.2450) . Finally, we find that the maximal expected revenue is ER(10) ≈ 1.7553 million dollars in mission interval [0,10] , and it is obtained when r = 1 and q ≈ 0.2449 , which means each mission repeated try succeeds with probability p = 1 − q ≈ 0.7551 .

7. Conclusions In this paper, we consider time redundant systems whose mission must be completed in a given mission window which is long enough for several mission repeated tries. In the main model, short mission repeated tries are assumed to fail definitely, and long ones are assumed to fail with decreasing probabilities; besides, a failed mission repeated try is followed by a repair duration before a new mission repeated try until the mission succeeds or stops. Models 1-3 are three special cases of the main model, Model 1 assumes failure rates of long mission repeated tries in the main model to be constant, Model 2 assumes repair durations in Model 1 to be zero, and Model 3 assumes that there is no short mission repeated try in Model 2. The main model and its three special cases are investigated by the Laplace transform method and 30

their relationships are also presented. To compare the Laplace transform method with the Monte Carlo simulation method and to consider different failure modes in our model, some extended discussions of the main model are presented. The optimal analysis for cost and revenue are also presented for the main model. Some numerical examples are provided to demonstrate our presented results and to illustrate the optimal analysis. Models in this paper can be applied to many practical systems in the field of reliability and some optimization problems considering different kinds of redundancy simultaneously can be studied in the future.

References Balakrishnan, N., Beutner, E., & Kamps, U. (2011). Modeling parameters of a load-sharing system through link functions in sequential order statistics models and associated inference. IEEE Transaction on Reliability, 60, 605-730. Barlow, R. E., & Proschan, F. (1996) Mathematical theory of reliability. Philadelphia: SIAM. Chambari, A., Rahmati, S. H. A., & Najafi, A. A. (2012). A bi-objective model to optimize reliability and cost of system with a choice of redundancy strategies. Computers & Industrial Engineering, 63, 109-119. Chang, G. J., Cui, L. R., & Hwang, F. K. (2000). Reliabilities for consecutive-k systems. Dordrecht: Kluwer Academic Publishers. Cherkesov, G. N. (1974). Reliability of technical systems with time redundancy. Moscow: Sov. Radio. Cohen, A. M. (2007). Numerical methods for Laplace transform inversion. New York: Springer. Cui, L. R., Lin, C., & Du, S. J. (2015). m-consecutive-k, l-out-of-n systems. IEEE Transaction on Reliability, 64, 386-393. Eryilmaz, S. (2017). The effectiveness of adding cold standby redundancy to a coherent system at system and component levels. Reliability Engineering & System Safety, 165, 331335. Eryilmaz S, Mahmoud B (2012). Linear m-consecutive-k, l-out-of-n: F system. IEEE Transactions on reliability, 61: 787-791. 31

Heidtmann, K. (1982). Time redundancy with applications to electronical devices. In Proc.5th Europ. Conf. Electrotechnics, EUROCON (vol. 82, pp. 188–191). Ivanov, D. V. (2012). Advanced sensors for multifunctional applications. Chemical Business, 26, 36-41. Jia, H. P., Ding, Y., Peng, R., & Song, Y. H. (2017). Reliability Evaluation for DemandBased Warm Standby Systems Considering Degradation Process. IEEE Transaction on Reliability, 66, 795-805. Kuhfittig, P. K. (2013). Introduction to the Laplace transform. New York: Springer. Kuznetsov, A. (2013). On the convergence of the Gaver--Stehfest algorithm. SIAM Journal on Numerical Analysis, 51, 2984-2998, 2013. Kvam, P. H., & Pena, E. A. (2006). Estimating load-sharing properties in a dynamic reliability system. Journal of the American Statistical Association, 100, 262-272. Langer, Y. (2001). System with time-redundancy. In Reliability and Maintainability Symposium, 2001. Proceedings. Annual. IEEE. (pp. 189-192). Levitin, G., Xing, L., & Dai, Y. (2013). Optimal sequencing of warm standby elements. Computers & Industrial Engineering, 65, 570-576. Lin, C., Cui, L. R., Coit, D. W., & Lv, M. (2016). Reliability modeling on consecutive- kr out-of- nr : F linear zigzag structure and circular polygon structure. IEEE Transaction on Reliability, 65, 1509-1521. Lin, Y. K., & Nguyen, T. P. (2018). Reliability evaluation of a multistate flight network under time and stopover constraints. Computers & Industrial Engineering, 115, 620-630. Linton, D. G., & Saw, J. G. (1974). Reliability analysis of the k-out-of-n: F system. IEEE Transaction on Reliability, 23, 97-103. Lisnianski, A., & Jeager, A. (2000). Time-redundant system reliability under randomly constrained time resources. Reliability Engineering & System Safety, 70, 157-166. Lisnianski, A., Levitin, G., & Ben-Haim, H. (2000). Structure optimization of multi-state system with time redundancy. Reliability Engineering & System Safety, 67, 103-112. Liu, H. M. (1998). Reliability of a load-sharing k-out-of-n:G system: Non-iid components with arbitrary distributions. IEEE Transaction on Reliability, 47, 279-284. 32

Nagel, J., Buchter, A., Xue, F., Kieler, O. F., Weimann, T., Kohlmann, J., & Berberich, P. (2013). Nanoscale multifunctional sensor formed by a Ni nanotube and a scanning Nb nanoSQUID. Physical Review B, 88, 064425. Pallaro, N., Visintainer, F., Darin, M., & Mosca, E. (2004). Multifunctional sensor to detect environmental parameters. In J. Valldorf, & W. Gessner (Eds), Advanced Microsystems for Automotive Applications (pp. 247-265). Berlin: Springer. Qiu, Q., Cui, L., & Gao, H. (2017a). Availability and maintenance modelling for systems subject to multiple failure modes. Computers & Industrial Engineering, 108, 192-198. Qiu, Q., Cui, L., Shen, J., & Yang, L. (2017b). Optimal maintenance policy considering maintenance errors for systems operating under performance-based contracts. Computers & Industrial Engineering, 112, 147-155. Rao, T. S., & Gupta, U. C. (2000). Performance modelling of the M/G/1 machine repairman problem with cold-, warm-and hot-standbys. Computers & Industrial Engineering, 38, 251267. Roozeboom, C. L., Hopcroft, M. A., Smith, W. S., Sim, J. Y., Wickeraad, D. A., Hartwell, P. G., & Pruitt, B. L. (2013). Integrated multifunctional environmental sensors. Journal of Microelectromechanical Systems, 22, 779-793. Shao, J. J., & Lamberson, L. R. (1991). Modeling a shared-load k-out-of-n:G system. IEEE Transaction on Reliability, 40, 205-209. Stehfest, H. (1970a). Algorithm 368: Numerical inversion of Laplace transforms. Communications of the ACM, 13, 47-49. Stehfest, H. (1970b). Remark on algorithm 368: Numerical inversion of Laplace transforms. Communications of the ACM, 13, 624-624. Ushakov, I. A. (1994). Handbook of reliability engineering. New York: Wiley. Wu, X. Y., & Hillston, J. (2015). Mission reliability of semi-Markov systems under generalized operational time requirements. Reliability Engineering & System Safety, 140, 122-129. Xu, Y. W., & H. T. Liao (2016). Reliability analysis and redundancy allocation for a one-shot system containing multifunctional components. IEEE Transaction on Reliability, 65, 10451057. 33

Yu, H., Yang, J., Lin, J., & Zhao, Y. (2017). Reliability evaluation of non-repairable phasedmission common bus systems with common cause failures. Computers & Industrial Engineering, 111, 445-457. Zhao, P., Zhang, Y. Y., & Chen, J. B. (2017). Optimal allocation policy of one redundancy in a n-component series system. European Journal of Operational Research, 257, 656-668. Zhao, P., Zhang, Y. Y., & Li, L. (2015). Redundancy allocation at component level versus system level. European Journal of Operational Research, 241, 402-411. Zhu, X. Y., Boushaba, M., Coit, D. W., & Benyahia, A. (2017). Reliability and importance measures for m-consecutive-k, l-out-of-n system with non-homogeneous Markov-dependent components. Reliability Engineering & System Safety, 167, 1-9.

34

Highlights: (1) (2) (3) (4)

Three time redundancy models are developed for a given mission window. The reliability analysis including optimal solutions is given under three models. The relationship among three models is considered. An algorithm for Monte-Carlo simulation is presented.

35

All Figures in the paper

mission window 2 1

α

α

observation station

horizontal line

orbit of the spaceflight

the earth

Figure 1. Mission window for data transmission of a TT&C system

N (t ) n n −1

a

stop

rq

q

1 1− r kn q

2 1 0

X1

Y1

Y2

X2

X3 mission window

τ

Xn

t

Figure 2. Mission success case for the main model.

N (t ) n n −1

2 1 0

a

X1

1

q

Y1

X2

stop

rq

r kn q

Y2

X3 mission window

Xn

Figure 3. Mission failure case for the main model.

36

X n +1

τ

t

N1 (t ) n n −1

stop

a

q

1 p

2 1 0

X1

Y2 X2 mission window

Y1

t

τ

Xn

Figure 4. Mission success case for sub-model 1.

N 2 (t ) n n −1

stop

a

q

1 p

2 1 0

X1

X2

Xn mission window

τ

t

τ

t

Figure 5. Mission success case for sub-model 2.

N 3 (t ) n n −1

stop

q

q p

2 1 0

X1

X2

Xn mission window

Figure 6. Mission success case for sub-model 3.

37

Figure 7. E[ N (t )] (Left) and H (t ) (Right) for the main model and its sub-models.

Figure 8. E[ N (t )] (Left) and H (t ) (Right) for the main model and its sub-models by the inverse Laplace method (solid lines) and by the Monte Carlo method (dotted lines).

38

Figure 9. The expected net revenue in interval [0,10] for different parameters q, r .

39

Table 1. Summary of the main model and its sub-models.

Failure probability of the

ith

Whether a repair duration exists

Models mission repeated try For X i ≤ a :

after a failed mission repeated try

1

Main model For X i > a : qki = r q

Yes, the durations Yi ~ G ( y )

ki

For X i ≤ a :

1 Yes, the durations Yi ~ G ( y )

Sub-model 1 For X i > a :

q

For X i ≤ a :

1

Sub-model 2

No For X i > a :

Sub-model 3

q

q

No

40

Acknowledgements

This work is supported by the National Natural Science Foundation of China under grant 71631001 and by China Scholarship Council. We also thank the associate editor and two referees’ valuable comments on improving the paper.

41