Reliability modeling with condition-based maintenance for binary-state deteriorating systems considering zoned shock effects

Accepted Manuscript Reliability modeling with condition-based maintenance for binary-state deteriorating systems considering zoned shock effects Guanzhou Wei, Xiujie Zhao, Shuguang He, Zhen He PII: DOI: Reference:

S0360-8352(19)30119-6 https://doi.org/10.1016/j.cie.2019.02.034 CAIE 5720

To appear in:

Computers & Industrial Engineering

Received Date: Revised Date: Accepted Date:

24 July 2018 7 January 2019 19 February 2019

Please cite this article as: Wei, G., Zhao, X., He, S., He, Z., Reliability modeling with condition-based maintenance for binary-state deteriorating systems considering zoned shock effects, Computers & Industrial Engineering (2019), doi: https://doi.org/10.1016/j.cie.2019.02.034

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.

Reliability modeling with condition-based maintenance for binary-state deteriorating systems considering zoned shock effects

Authors Guanzhou Weia (first author); Email: guanzhou [email protected] Xiujie Zhaob (coauthor); Email: [email protected] Shuguang Hea∗ (corresponding author); Email: [email protected] Zhen Hea (coauthor); Email: [email protected] a

College of Management & Economics, Tianjin University, Tianjin, China

Postal address: Quality Management and Quality Engineering Lab, 25 Building, Tianjin University, 92 Weijin Road, Nankai District, Tianjin, China b

Department of Systems Engineering and Engineering Management, City University of Hong

Kong, Kowloon, Hong Kong Postal address: P7306, Quality and Reliability Engineering Lab, Yeung Kin Man Academic Building (AC1), City University of Hong Kong, 83 Tat Chee Avenue, Kowloon Tong, Hong Kong.

Declarations of interest: none

1

Reliability modeling with condition-based maintenance for binary-state deteriorating systems considering zoned shock effects

Abstract Due to the necessity of usage, industrial systems often work and deteriorate in a dynamic environment, where some deleterious factors can be regarded as external shocks. In this study, we consider a binary-state deteriorating system under zoned shock effects. System internal degradation is modeled by a two-phase Wiener process, which has a larger shift and diffusion parameter when the system transfers from the normal state to the weakened state. Shock effects are comprehensively incorporated in the deteriorating system, where the imperfect shocks are added to the zoned shock effects. Further, a new condition-based maintenance (CBM) strategy is built to improve the system performances. An optimal action (i.e., no action, imperfect repair, preventive replacement or corrective replacement) is determined based on the system state by minimizing the average long-run cost rate. Finally, a numerical example is presented to illustrate the applicability of the proposed reliability model with CBM.

Keywords: Two-phase degradation path; Zoned shock effects; Condition-based maintenance

1

1

Introduction Complex mechanical devices and engineering infrastructures (e.g., micro-electro-mechanical

systems (MEMS), oil pipelines, dams etc.) generally deteriorate with functioning and aging increasingly. System deterioration results from their underlying degradation mechanisms like wear, fatigue, corrosion and so on, which can ultimately cause systems to malfunction. The functional degradation-based failure is generally called the soft failure (Ye et al., 2012). Commonly, such deteriorating systems operate in a dynamic environment to fulfill missions and some unexpected events may occur to the system from time to time. For example, a river-crossing bridge is designed for transportation, withholding fluctuation of workloads, and it is subject to some extreme events like floods and hurricanes. When a system suffers from a sudden and unexpected accident such as extreme load, stress, and voltage, it is most likely to be destroyed. Such failures caused by traumatic events or shocks is often referred to as the hard failure. In most cases, soft failures and hard failures are competing, which is known as degradation-threshold-shock (DTS) model (Lehmann, 2009; Huynh et al., 2012). In DTS models, soft failures occur when the system degradation level exceeds a predetermined threshold. Due to the internal randomness among individuals and the variability from ambient environment factors like temperature, humidity etc., system degradation is often modeled by stochastic processes. Commonly used stochastic processes include Gamma process, inverse Gaussian process and Wiener process (Lawless and Crowder, 2004; Wang and Xu, 2010; Peng, 2015; Huang et al., 2015). Among them, Wiener process is extensively employed in degradation analysis, see Zhang et al. (2018) for an overview. It has satisfactory mathematical tractability and avoids undesirable model complexity, providing an effective way to evaluate the reliability of deteriorating systems (Zhang et al., 2017) and support further decision-making on maintenance (Liu et al., 2016). As for hard failures caused by external shocks, they often occur to systems when a single shock hits a certain level. Such shocks with a magnitude higher than the threshold could ruin a system directly. Whereas not every shock has the ability to devastate systems catastrophically. Recently, zoned shock effects have been considered by dividing random shocks into three zones based on their magnitudes, i.e., the safety zone, the damage zone, and the fatal zone. Peng et al. (2010) described random shocks with a smaller magnitude as not destroying systems catastrophically but bringing instantaneous degradation increments to accelerate the soft failure. Jiang et al. (2015) found that micro-engines had the ability to resist small shocks owing to the material strength. Apart from this, Nakagawa (2007) argued that there should exist an interme-

2

diate zone for imperfect shocks between the safety zone and the damage zone, where imperfect shocks only have a certain probability to damage the system. However, there is no research that considers the imperfect shocks in reliability modeling. Thus, we attempt to cover the effect of imperfect shocks in this study. While most existing studies on the DTS model assumed that systems had a consistent deterioration characteristic, systems with a two-phase deterioration path commonly exist in industrial practices. Bae and Kvam (2006) found that the degradation path of Plasma Display Panels (PDPs) exhibited a two-phase pattern during burn-in tests. The light intensity decreased rapidly for a short period and then followed its inherent path with a relatively slow degradation rate. Another example is the fatigue crack propagation of mental components. The propagation rate of cracks increases sharply when the crack length reaches a critical level (Zhang et al., 2016). Despite the existence of these realistic motivations, only limited researches focus on this issue. Yang et al. (2017a,b) considered two-stage deterioration systems in which external shocks can only devastate system directly. Zhang et al. (2016) adopted the Wald’s equation to approximately model the overall degradation of a two-state system, which integrates degradation caused by damage zone shocks instantaneously. Nevertheless, safety zone shock effects are not discussed and no researches comprehensively incorporate zoned shock effects to binary-state deteriorating systems owing to the challenging complexity of modeling. Additionally, if zoned shock effects could be investigated in binary-state systems thoroughly, one could have a better insight into system failures so as to make more efficient maintenance decisions. Motivated by this need, we discuss zoned shock effects in binary-state deteriorating systems comprehensively. With the advancement of technologies, degrading systems are nowadays equipped with monitoring mechanisms such as sensors, meters, and computational devices, making it possible to assess system health conditions online. For this reason, condition based maintenance (CBM) policies for such deteriorating systems have received intensive attention over the past years. One can see Alaswad and Xiang (2017) and Keizer et al. (2017) for overviews. CBM considers the real-time condition of systems such as the system state and degradation levels through periodic inspections or continuous monitoring. When failures are detected or the degradation level exceeds the predetermined preventive maintenance (PM) threshold, corrective maintenance (CM) or PM are employed respectively to improve system performances. Castro et al. (2015) proposed a CBM strategy with periodic inspections for degrading systems under random shocks. A preventive replacement was performed when the degradation exceeded the predetermined threshold and a corrective replacement was implemented when the system was failed. Yang et al. 3

(2017b) considered different preventive replacement thresholds in two-state deteriorating systems combining corrective replacement and preventive replacement in CBM. Nevertheless, few studies modeled imperfect PM or imperfect repair in CBM for aforementioned binary-state deteriorating systems since imperfect repair aggravates the complexity of maintenance modeling. While Do and B´erenguer (2012) and Do et al. (2015) investigated imperfect PM combined with perfect PM in CBM for deteriorating systems through simulation methods, researches that model CBM by integrating imperfect repair is relatively lacking, especially for binary-state deteriorating systems with zoned shock effects. The modeling is very challenging and this study intends to address this issue. Distinct from existing works, we address three main challenging problems in this study. 1) We add imperfect shocks in the reliability modeling, which extends zoned shock effects. 2) Zoned shock effects are comprehensively discussed in the binary-state system including safety zone shocks, imperfect shocks, damage zone shocks, and fatal shocks. 3) A new CBM strategy is proposed by combining imperfect repair, preventive replacement and corrective replacement together for the binary-state system under zoned shock effects. The remainder of this paper is organized as follows. Section 2 presents the definition of the system state and depicts the arrival of external shocks. Section 3 discusses the system reliability modeling including system internal degradation and zoned shock effects. In Section 4, we propose a new CBM strategy for the binary-state system, and the optimization of CBM is discussed. A numerical example is presented in Section 5 to illustrate the applicability of the reliability model with the optimal CBM. Finally, we give concluding remarks in Section 6.

2

System description In this section, a binary-state deteriorating system under random shocks is considered. We

define two states of the system and model the arrivals of random shocks to the system under its two states respectively.

2.1

Binary-state system

For deteriorating mechanical systems, their degradation level increases gradually and the structural strength becomes weaker, reducing the resistance ability of systems to external shocks. The overall degradation level of a system can reflect its health condition to some extent. Based on the overall degradation level, we define a functioning system with two states, i.e., the normal 4

state and the weakened state. As shown in Fig. 1, the employable system is in the weakened state when its overall degradation level X(t) is between the predetermined weakened state threshold H1 and the degradation-based failure threshold H2 . Similarly, the normal state is when the overall degradation level is smaller than H1 . We use the random variable Tw to denote the duration transferring from the normal state to the weakened state. Specifically, Tw is the first passage time when the overall degradation level reaches H1 given X(0) = 0 and tw is its realization. Note that in Fig. 1, shock effects on system deterioration are shown in advance, which will be discussed in Section 3.2.

X (t) H2

Normal state

H1

Weakened state

tw

t

Fig. 1. Binary-state system and its degradation path

2.2

Arrival of external shock process

We adopt an homogeneous Poisson process (HPP) {N (t), t ≥ 0} with rate λ to model random shocks, where N (t) denotes the number of all shocks in the time interval (0, t]. For any t ≥ 0, N (t) is a non-negative random variable following a Poisson distribution with the expression of: P(N (t) = n) =

e−λt (λt)n , n = 0, 1, · · · . n!

(1)

Let Wi be the magnitude of the ith shock si and we assume Wi s(i = 1, 2, · · · ) are identically independent random variables with a cumulative distribution function (CDF) FW (w), which implies random shocks are independent in the magnitude dimension. As shown in Fig. 2, random shocks are classified into four zones based on their magnitude W , i.e., the fatal zone (FZ), the damage zone (DZ), the probability damage zone (PDZ), and the safety zone (SZ). In general, when a system transfers from the normal state to the weakened state, its resistance to external shocks declines. A reasonable assumption is that a system in the weakened state can only resist relatively small shocks. In our proposed model, it can be represented as a reduction of the shock 5

zoned thresholds shown in Fig. 2, where D30 = D3 − ∆, D20 = D2 − ∆, D10 = D1 − ∆ and ∆ is a constant depending on the characteristics of a system such as the material type and the design lifespan. X (t) H2

Y2( 2) Y2( 3) H1 Y1( 2) Y1( 3)

Y1( 3)

t1 t2

t3

fatal zone

t1

t2

t3

t4

t5

tw

t6

t7

probability damage zone

damage zone

t1

t

safety zone

W

D3

∆

D2

∆

D1

∆

t1 t2

t3

t1

t2

t3

t4

t5

tw

t6

t7

t1

t

Fig. 2. Zoned shocks under the binary-state system

Table 1. Classification of shock zones in the binary-state system (a) Classification of shock zones in the normal state system Shock zone FZ DZ PDZ SZ

Magnitude [D3 ,+∞) [D2 ,D3 ) [D1 ,D2 ) [0,D1 )

Probability of random shocks belonging to the zone 1 − FW (D3 ) FW (D3 ) − FW (D2 ) FW (D2 ) − FW (D1 ) FW (D1 ) − FW (0)

(b) Classification of shock zones in the weakened state system Shock zone FZ DZ PDZ SZ

Magnitude Probability of random shocks belonging to the zone [D3 − ∆,+∞) 1 − FW (D3 − ∆) [D2 − ∆,D3 − ∆) FW (D3 − ∆) − FW (D2 − ∆) [D1 − ∆,D2 − ∆) FW (D2 − ∆) − FW (D1 ) − ∆ [0,D1 − ∆) FW (D1 − ∆) − FW (0)

The probabilities of the ith shock si falling into the SZ, the PDZ, the DZ, and the FZ are p1 , p2 , p3 and p4 expressed in Table 1 (a). The counting numbers of shocks falling into the four zones in the time interval (0, t] are denoted by N1 (t), N2 (t), N3 (t) and N4 (t) respectively. Since p1 + p2 + p3 + p4 = 1, {Ni (t), t ≥ 0, i = 1, 2, 3, 4} are independent HPPs with rate λpi according to the decomposition theorem of the HPP (Ross, 2014). Thus, the probability mass function (PMF) of Ni (t) can be obtained as: P{Ni (t) = n} =

e−λpi t (λpi t)n , n = 0, 1, · · · , i = 1, 2, 3, 4. n! 6

(2)

Similarly, when the system is in the weakened state, the probabilities of the ith shock si falling into the changed shock zones (i.e., the SZ, the PDZ, the DZ and the FZ under the weakened state system) are p01 , p02 , p03 and p04 expressed in Table 1 (b) respectively. Let {Ni0 (t), t ≥ tw , i = 1, 2, 3, 4} denote the processes of random shocks arriving at the four shock zones when the system is in the weakened state. The four processes follow HPPs with rate λp01 , λp02 , λp03 , and λp04 , and the PMFs of Ni0 (t) (i = 1, 2, 3, 4), denoting the total numbers of random shocks arriving at the four shock zones in the time interval [tw , t), is the expression of: 0

P{Ni0 (t)

e−λpi (t−tw ) [λp0i (t − tw )]n = n} = , n = 0, 1, · · · , i = 1, 2, 3, 4. n!

(3)

Remark 1. The Poisson process has the independent inter-arrival time between events following the exponential distribution, which is often used to describe shocks arriving in reliability modeling and customer arrivals in queuing theory. Such property is useful since the exponential distributed inter-arrival time provides large flexibility in realistic applications. When choosing appropriate parameters, the assumption could be applied to a variety of engineering and managerial problems.

3

Reliability modeling For the proposed binary-state system, both external shocks and system internal degradation

are related to the system state, where main assumptions are summarized as follows: 1) The internal degradation path can be modeled by Wiener process and its degradation rate in the weakened state is greater than that in the normal state. 2) The damage Y (2) caused by random shocks in the PDZ and Y (3) caused by random shocks in
the DZ are random variables following the normal distributions, i.e., Y (2) ∼ N µp , σp2 and Y (3) ∼ N (µd , σd2 ), where µp and µd are the expected additional increments on system internal degradation caused by each shock in the DZ and the PDZ respectively. 3) The probability of a shock in the PDZ leading damage is in proportion to the shock magnitude. A larger magnitude of the imperfect shock will lead to a larger probability of damage effects to the system. Apart from this, if a random shock with a magnitude w in the PDZ and the DZ causes an additional degradation y, we measure it as y = α1 w when a system is in the normal state and α2 w in the weakened state. Remark 2. After suffering a series of random shocks, a system becomes more sensitive to degradation and shifts to an accelerating degradation period with a faster degradation rate (Rafiee 7

et al., 2014). Additionally, in order to proximately provide a monotone increasing degradation path, the shift parameter is far greater than the diffusion parameter in the Wiener process, i.e, µp  σp and µd  σd . Such approximation is valid due to the small probability of negative degradation increments. Remark 3. Commonly, the damage caused by random shocks is complicated. It results from not only various characteristics of random shocks (e.g. size, inter-arriving time, and hitting direction) but also the system conditions including surface roughness, fatigue level, etc. Considering all the contributions from these factors, random shock effects can be approximately modeled by the normal distribution. Remark 4. The probability of a shock leading to system damage is zero in the SZ while one in the DZ. One reasonable assumption is that it increases as shock magnitudes rise. Coefficients α1 , α2 are constants and they may be experimentally determined in a specific application or based on large amounts of historical knowledge. Note that α2 > α1 implies systems in the weakened state are more vulnerable to random shocks.

3.1

System internal degradation

System internal degradation over time can be modeled by Wiener process {Xc (t), t ≥ 0} with a drift parameter µc and a diffusion parameter σc > 0, where Xc (t) = µc t + σc B(t) and B(t) is a standard Brownian motion. It satisfies: 1) Xc (0) = 0 2) {Xc (t), t ≥ 0} has stationary and independent increments 3) {Xc (t), t ≥ 0} is normally distributed with mean µc t and variance σc2 t for any t > 0. As mentioned in model assumptions, the system internal degradation corresponds to the system state. It can be expressed as: ( Xc (t) =

Xc1 (t) = µc1 t + σc1 B(t),

0 ≤ t < tw

Xc2 (t) = Xc1 (tw ) + µc2 (t − tw ) + σc2 B(t − tw ),

t ≥ tw > 0,

(4)

where µc2 > µc1 > 0 and they can be viewed as mean degradation rates of the system in two states respectively. The diffusion parameter σc also changes with the system state, i.e., σc2 > σc1 > 0, which implies that systems in the weakened state have larger volatility in deterioration. From Eq. 2 (4), for any t2 < t1 < tw , Xc1 (t2 ) − Xc1 (t1 ) ∼ N (µc1 (t2 − t1 ), σc1 (t2 − t1 )) and Xc2 (t2 ) − Xc2 (t1 ) ∼ 2 N (µc2 (t2 − t1 ), σc2 (t2 − t1 )) for t2 > t1 > tw .

8

3.2

Modeling zoned shock effects

Random shocks may either ruin system structures directly or bring some damage to them gradually, where damage can be described by additional increments to system internal deterioration. Shocks in the DZ contribute to soft failures by adding abrupt increments to system internal deterioration, and imperfect shocks have a probability p(w) to bring an additional degradation level to accelerate system soft failures. FZ shocks can devastate the system structures catastrophically, while shocks in the SZ are considered harmless owing to the system strength. As mentioned earlier, fatal shocks come to a system according to an HPP with rate λp4 when the system is in the normal state and rate λp04 in the weakened state. We denote Th as the hard failure time caused by fatal shocks. When a system is in the normal state, the survival function of Th conditional on a given tw using Eq. (2) is obtained as: F¯Th |tw >t (t) = P(Th > t|tw > t) = P{N4 (t) = 0} = e−λp4 t .

(5)

For a system in the weakened state, the survival function of Th conditional on a given tw with Eqs. (2) and (3) is derived as: F¯Th |tw ≤t (t) = P(Th > t|tw ≤ t) = P{N4 (tw ) = 0, N40 (t) = 0} 0

0

0

= e−λp4 tw · e−λp4 (t−tw ) = eλ[(p4 −p4 )tw −p4 t] .

(6)

Since each shock in the DZ brings an instantaneous degradation shift, the cumulative damage caused by random shocks in the DZ can be modeled by a two-stage compound Poisson process (3)

{Xs (t), t ≥ 0} shown as follows:

Xs(3) (t) =

 N3 (t) X (3)   (3)  X (t) = Yi ,    s1

0 ≤ t < tw

i=0

N30 (t)  X (3)   (3) (3)  X (t) = X (t ) + Yj ,  w s2 s1 

t ≥ tw > 0,

(7)

j=0

(3)

where Yi

is an additional increment caused by the ith shock in the DZ when a system is in the (3)

normal state and Yi

(3)

∼ N (µd1 , σd2 ). Yj

is an additional increment caused by the jth shock in (3)

the DZ when a system is in the weakened state and Yj

9

∼ N (µd2 , σd2 ). Further, the expected

damage µd1 is obtained as: Z

D3

α1 wdFW |W ∈(D2 ,D3 ) (w) = α1 E [W |W ∈ (D2 , D3 )] ,

µd1 =

(8)

D2

where FW |W ∈(D2 ,D3 ) (w) is the conditional CDF of the shock magnitude in the DZ. Similarly, µd2 = α2 E[W |W ∈ (D20 , D30 )] For a shock in the PDZ, we assume that the probability of shocks leading to additional degradation increases linearly as shock magnitudes rise and it is defined as: p(w) =

w − D1 , w ∈ (D1 , D2 ). D2 − D1

(9)

This definition is consistent with Finkelstein and Zarudnij (2001), which implies that the probability of shocks causing damage increases as their magnitudes approach to D1 . We measure the expected damage caused by each imperfect shock for a system in the normal state by: Z

D2

α1 w · p(w)dFW |W ∈(D1 ,D2 ) (w)

µp1 = D1

=

α1 {V[W |W ∈ (D1 , D2 )] + {E[W |W ∈ (D1 , D2 )]}2 } − α1 D1 E[W |W ∈ (D1 , D2 )] . D2 − D1

(10)

Likewise, the expression of µp2 can be derived. Furthermore, similar to the DZ shocks, the cumulative shock damage caused by random shocks in the PDZ is given by:

Xs(2) (t) =

 N2 (t) X (2)   (2)  X (t) = Yi ,    s1

0 ≤ t < tw

i=0

N20 (t)  X (2)   (2) (2)  X (t) = X (t ) + Yj ,  w s2 s1 

t ≥ tw > 0,

(11)

j=0

(2)

where Yi

is the additional degradation caused by the ith imperfect shock when a system is in
(2) (2) the normal state and Yi ∼ N µp1 , σp2 . Yj is the additional degradation caused by the jth
(2) imperfect shock when a system is in the weakened state and Yj ∼ N µp2 , σp2 .

3.3

System reliability function

The overall degradation X(t) is the sum of the system internal degradation Xc (t), the addi(2)

(3)

(2)

(3)

tional damage Xs (t) and Xs (t), i.e., X(t) = Xc (t)+Xs (t)+Xs (t). The soft failure occurs to a system once the overall degradation level exceeds the predetermined degradation-based thresh10

old H2 . Note that it can only occur when a system is in the weakened state. We use Td to denote the soft failure time, and the survival function of Td conditional on a given tw is the expression of: F¯Td |tw ≤t (t) =

∞ X ∞ X

P{Td > t|N20 (t) = k, N30 (t) = l} · P{N20 (t) = k} · P{N30 (t) = l},

(12)

k=0 l=0

where P{N20 (t) = k}, P{N30 (t) = l} can be obtained by Eq. (3), and P{Td > t|N20 (t) = k, N30 (t) = l} is derived as: P{Td > t|N20 (t) = k, N30 (t) = l} = P{X(t) − X(tw ) < H2 − H1 |N20 (t) = k, N30 (t) = l}   k l   X X (2) (3) = P Xc (t) − Xc (tw ) + Yj + Yj < H2 − H1   j=0 j=0   H2 − H1 − µc2 (t − tw ) − kµp2 − lµd2  q = Φ . 2 (t − t ) + k 2 σ 2 + l2 σ 2 σc2 w p d

(13)

If a system has survived by time t, it is ensured that both the soft failure and the hard failure did not happen before time t. Based on the total probability theorem and Eqs. (5),(6) and (12), the system reliability is obtained as: Z R(t) =

∞

F¯Th |u>t (t)dFTw (u) +

t

t

Z

F¯Th |u≤t (t; u) · F¯Td |u≤t (t; u)dFTw (u),

(14)

0

where FTw (·) is the CDF of Tw . See the derivation of FTw (·) in Appendix A.

4

CBM strategy In this study, we propose a new CBM strategy with periodic inspections. As shown in Fig. 3,

systems are inspected at time epochs kτ (k = 1, 2, · · · ), and maintenance decisions are based on the system state at each inspection. Main assumptions of the proposed CBM are as follows: 1) The system is repairable and the repair action can be performed only once between two consecutive replacements. 2) Each inspection is instantaneous and nondestructive. Repair and replacement time are negligible compared with the system lifespan. The system is not self-announcing, and the failed system remains idle until the next inspection. 3) At the ith inspection, (i) if the system still works and is in the normal state, we leave it out 11

and take no maintenance action. (ii) if the system is detected to fail either due to soft failures or hard failures, corrective replacement is conducted. (iii) if the system is in the weakened state and its overall degradation level is smaller than the preventive replacement threshold H, we carry out an imperfect repair action to restore the degradation level to H1 , where H1 < H. (iv) if the system is still functioning and its overall degradation level is between H and H2 , preventive replacement is performed, where H2 > H. X (t) H2

no action

H

imperfect repair

Y4( 3)

Y3( 3)

preventive replacement

H1

Y2( 3)

corrective replacement

Y1( 2) Y1( 3)

Y1( 3)

τ

4τ

3τ

2τ

7τ

6τ

5τ

t

W

fatal zone

∆

D3

damage zone

∆

D2

probability damage zone

D1

∆

safety zone

t1

t2

t3

t4

t5

t6

t7

t1

t2

t

Fig. 3. The proposed CBM strategy

Remark 5. Practically, the repair action can only be performed for limited times since repairing is harmful to functioning systems. In our proposed CBM, if repair actions are not restricted, short inspection periods lead to infinite times of the imperfect repair. Apart from this, due to the challenging complexity of modeling, we adopt once repair limitation in Rafiee et al. (2015) to simplify our further derivation. Remark 6. Each inspection incurs an inspection cost CI . If the system is implemented with corrective replacement, it costs more than the preventive replacement in general. We assume the corrective replacement cost CR is larger than the preventive replacement cost CP . Further, the imperfect repair cost CIR is smaller than the preventive replacement cost CP while larger than the inspection cost CI , i.e., CR > CP > CIR > CI . In the proposed CBM strategy, a renewal cycle is defined as the whole lifetime of a functioning system. It is a period either from the system installation to the first replacement or two consecutive replacement actions, where the replacement refers to as preventive replacement or corrective 12

replacement. Each cycle incurs a related cost, and these costs along with the cycles form a renewal process. To evaluate the CBM performance, the average long-run total maintenance cost rate CR is adopted, and two decision variables corresponding to CR are the inspection period τ and the preventive replacement threshold H. Therefore, we have: C(t) E(T C) = , t→∞ t E(P )

CR(τ, H) = lim

(15)

where C(t) is the total maintenance cost by time t. E(T C) is the expected total maintenance cost of a renewal cycle and E(P ) is the expected cycle length. In order to evaluate CR, we derive the probabilities of performing different types of maintenance actions firstly.

4.1

Probability of performing imperfect repair

The imperfect repair action takes place when a weakened state system still works and its overall degradation level is smaller than the preventive replacement threshold H. As shown in Fig. 4, the repair is performed at the ith inspection, and it means that the system is in the normal state at the (i − 1)th inspection while in the weakened state at the ith inspection, i.e., the transition time Tw is between (i − 1)τ and iτ . We denote NIR as the counting number of inspections until performing the imperfect repair. In order to simplify the following derivation, we define FX (x; t1 , t2 ) as the CDF of the overall degradation changes in the time interval (t1 , t2 ) under a weakened state system expressed as: FX (x; t1 , t2 ) = P{X(t2 ) − X(t1 ) ≤ x} =

∞ X ∞ X

P{N20 (t2 ) − N20 (t1 ) = k} · P{N30 (t2 ) − N30 (t1 ) = l}

k=0 l=0

·P{X(t2 ) − X(t1 ) ≤ x|N20 (t2 ) − N20 (t1 ) = k, N30 (t2 ) − N30 (t1 ) = l} 0 0 ∞ X ∞ X e−λp2 (t2 −t1 ) [λp02 (t2 − t1 )]k e−λp3 (t2 −t1 ) [λp03 (t2 − t1 )]l = · k! l! k=0 l=0 ( ) k l X X (2) (3) ·P Xc2 (t2 ) − Xc2 (t1 ) + Yj + Yj < x

j=0

=

∞ X ∞ X

−λp02 (t2 −t1 )

e

[λp02 (t2

k

− t1 )]

k!

k=0 l=0



·

e

j=0 −λp03 (t2 −t1 )

[λp03 (t2 − t1 )]l l! 

x − µc2 (t2 − t1 ) − kµp2 − lµd2  · Φ q , 2 2 2 2 2 σc2 (t2 − t1 ) + k σp + l σd

13

(16)

note that where t2 > t1 > tw and FX (x; t1 , t2 ) can only be usable when the system is in the weakened state. Then the PMF of NIR is derived as: Z

iτ

P {X(iτ ) − X(u) < H − H1 } · P{N4 (u) = 0} · P{N40 (iτ ) = 0}dFTw (u)

P(NIR = i) = (i−1)τ iτ

Z

0

FX (H − H1 ; iτ, u) · e−λp4 (iτ −u) dFTw (u).

=

(17)

(i−1)τ

X (t ) H2 H

H1

( i - 1)t

tw

it

t

Fig. 4. The imperfect repair at the ith inspection

4.2

Probability of performing preventive replacement

Preventive replacement is performed when the overall degradation level of a functioning system is between H and H2 , where H2 > H. Since the discussed system is repairable, it may undergo one imperfect repair action or no repair action before preventive replacement. Thus, the probability of performing preventive replacement at the ith inspection depends on whether the repair action has been done or not, and there are two scenarios: 1) Scenario 1: no imperfect repair action is conducted before preventive replacement. As shown in Fig. 5, at the (i − 1)th inspection, the system is in the normal state, while the weakened state is detected at the ith inspection. Meanwhile, its overall degradation level is between H and H2 . We use P(NP = i, NIR > i − 1) to denote the probability that the system is preventively replaced at the ith inspection without having undergone imperfect repair by the ith

14

inspection. With Eqs. (2), (3) and (16), the probability can be expressed by: iτ

Z P(NP = i, NIR > i − 1) =

P {H − H1 < X(iτ ) − X(u) < H2 − H1 } · P{N4 (u) = 0} (i−1)τ

· P{N40 (iτ ) = 0}dFTw (u) iτ

Z =

0

[FX (H2 − H1 ; iτ, u) − FX (H − H1 ; iτ, u)] · e−λp4 (iτ −u)

(i−1)τ

· e−λp4 u dFTw (u).

(18)

X (t )

X (t )

( i - 1)t

tH

tw

it

H2

H2

H

H

H1

H1

tw

t

Fig. 5. Scenario 1 for preventive replacement at the ith inspection

jt

( i - 1)t

tH

it

t

Fig. 6. Scenario 2 for preventive replacement at the ith inspection

2) Scenario 2: one imperfect repair action is conducted before preventive replacement. As shown in Fig. 6, since an imperfect repair action has been employed at the jth inspection, the overall degradation level of a repaired system begins with H1 at time jτ . The system is preventively replaced at the ith inspection when its overall degradation is between H and H2 , meanwhile, the degradation level is smaller than H at the (i − 1)th inspection in order to avoid preventive replacement happening at the (i − 1)th inspection. We denote P(NP = i|NIR = j) as the probability of performing preventive replacement at the ith inspection conditional on one imperfect repair action taking place at the jth inspection. It is derived as: Z

iτ

P(NP = i|NIR = j) =

P {X(iτ ) − X(u) < H2 − H} · P{N40 (iτ − jτ ) = 0}dFTH (u)

(i−1)τ −λp04 (iτ −jτ )

Z

iτ

FX (H2 − H; iτ, u)dFTH (u),

=e

(19)

(i−1)τ

where j ≤ i − 1 and TH is the first passage time when the degradation level reaches H given the initial degradation level H1 at time jτ . FTH (·) is the CDF of TH and we derive it in Appendix B.

15

Combining the two scenarios, we obtain the probability of performing preventive replacement at the ith inspection as:

P(NP = i) = P(NP = i, NIR > i − 1) +

i−1 X

P(NP = i|NIR = j) · P(NIR = j).

(20)

j=1

The analytical expression of Eq. (20) can be derived by Eqs. (17)-(19).

4.3

Probability of performing corrective replacement

The system undergoes corrective replacement when a failure is detected at the ith inspection caused by either fatal shocks or over-threshold deterioration. Similar to the preventive replacement situation, the probability of performing corrective replacement at the ith inspection depends on whether the repair action has been done or not. Therefore, there are two scenarios: 1) Scenario 1: no imperfect repair action is conducted before corrective replacement. As shown in Fig. 7, corrective replacement takes place at the ith inspection. The system is in the normal state at the (i − 1)th inspection, while it fails when fatal shocks occur between the (i − 1)th inspection and the ith inspection or the overall degradation level is higher than the degradation-based failure threshold H2 at the ith inspection. We use P(NR = i, NIR > i − 1) to denote the probability that corrective replacement is performed at the ith inspection, meanwhile no repair action has been employed by the ith inspection. Then we have: Z

iτ

P(NR = i, NIR > i − 1) =

{1 − P{X(iτ ) − X(u) < H2 − H1 }P{N4 [u − (i − 1)τ ] = 0} (i−1)τ

P{N40 (iτ ) = 0}} · P{N4 [(i − 1)τ )] = 0}dFTw (u) Z iτ 0 −λp4 (i−1)τ =e {1 − FX (H2 − H1 ; iτ, u) · e−λp4 [u−(i−1)τ ] · e−λp4 (iτ −u) } (i−1)τ

dFTw (u).

(21)

2) Scenario 2: one imperfect repair action is conducted before corrective replacement. As shown in Fig. 8, at the ith inspection, the system undergoes corrective replacement conditional on a repair action having been employed at the jth inspection. Thus, the repaired system deteriorates with the initial degradation level H1 at time jτ , and its overall degradation level ends with a value smaller than H at the (i − 1)th inspection as preventive replacement is not allowed in this scenario. Moreover, the system is still functioning at the (i − 1)th inspection, while a soft failure or a hard failure takes place between the (i − 1)th and the ith inspection. 16

We use P(NR = i|NIR = j) to denote the probability of performing corrective replacement at the ith inspection given that an imperfect repair action has been performed at the jth inspection (j ≤ i − 1), and we have: Z

iτ

{1 − P{X(iτ ) − X[(i − 1)τ ] < H2 − H1 }

P(NR = i|NIR = j) = (i−1)τ

P{N40 (iτ ) − N40 [(i − 1)τ ] = 0}} · P{N40 [(i − 1)τ ] − N40 (jτ ) = 0} dFTH (u) Z iτ

=e

−λp04 (i−j−1)τ

n o −λp04 τ 1 − FX (H2 − H1 ; iτ, (i − 1)τ )e dFTH (u), (22)

(i−1)τ

Combining the two scenarios, the probability of performing a corrective replacement action at the ith inspection can be obtained as:

P(NR = i) = P(NR = i, NIR > i − 1) +

i−1 X

P(NR = i|NIR = j) · P(NIR = j).

(23)

j=1

The analytical expression of Eq. (23) can be obtained by Eqs. (17), (21) and (22).

X (t)

X (t)

( i −1) τ

H2

H2

H

H

H1

H1

iτ

tw

( i −1) τ

t

iτ

tw

Fig. 7. Scenario 1 for corrective replacement at the ith inspection

X (t)

X (t)

( j −1) τ

tw

jτ

( i −1) τ

tH

iτ

H2

H2

H

H

H1

H1

t

( j −1) τ

tw

jτ

( i −1) τ

Fig. 8. Scenario 2 for corrective replacement at the ith inspection

17

iτ

4.4

Optimization of the CBM strategy based on CR

Since a renewal cycle ends with preventive replacement or corrective replacement, the expected renewal cycle length is derived as:

E(P ) = =

∞ X i=1 ∞ X

iτ · {1 − {1 − P(NP = i)} · {1 − P(NR = i)}} iτ · [P(NP = i) + P(NR = i)],

(24)

i=1

where the last step is valid because preventive replacement and corrective replacement are mutually exclusive, i.e., P(NP = i) · P(NR = i) = 0. Additionally, the analytical expression of E(P ) can be obtained by Eqs. (20) and (23). In a renewal cycle, total expected costs include three parts, i.e., E(T C) = E(IC) + E(RC) + E(IRC), where IC, RC, and IRC are the total inspection cost, the total replacement cost and the total imperfect repair cost in a renewal cycle respectively. Considering the expected inspection cost corresponds to the expected number of inspections in a renewal cycle, E(IC) is as follows: E(IC) = CI

∞ X

i[P(NP = i) + P(NR = i)].

(25)

i=1

The replacement cost is incurred by preventive replacement or corrective replacement in a renewal cycle. A renewal cycle expires at the ith inspection when a replacement action takes place, so the expected replacement cost is obtained as:

E(RC) = CP

∞ X

P(NP = i) + CR

i=1

∞ X

P(NR = i).

(26)

i=0

The imperfect repair takes place before a replacement action. If a system is performed a replacement action at the ith inspection, the imperfect repair cost is incurred when a repair action takes place before the ith inspection. Considering no repair action can be performed if replacement takes place at the first inspection, the imperfect repair cost is derived as:

E(IRC) = CIR = CIR

∞ X [P(NP = i, NIR ≤ i − 1) + P(NR = i, NIR ≤ i − 1)] i=2 ∞ X i−1 X

P[(NP = i|NIR = j) · P(NIR = j) + P(NR = i|NIR = j) · P(NIR = j)].

i=2 j=1

(27) 18

The analytical results of Eqs. (25)-(27) can be obtained by Eqs. (17), (19), (20), (22) and (23). As mentioned previously, the inspection period τ and the preventive replacement threshold H are two main factors corresponding to the CBM efficiency. Therefore, with the help of CR(τ, H) = E(T C) , E(P )

we can explore the optimal τ ∗ and H ∗ by solving the following programming problem. (τ ∗ , H ∗ ) = arg min CR(τ, H) τ,H

(

τ1 < τ < τ2

s.t. H1 < H < H2 ,

(28)

where τ1 and τ2 are constraints in a practical application. For the complexity of the subject function CR(τ, H), it is extremely difficult, if possible, to obtain the analytical solution of τ ∗ and H ∗ . Accordingly, we illustrate the optimal solution through a numerical example in Section 5.

5

Numerical example - a MEMS application To illustrate the applicability of our proposed reliability model with CBM, we take a MEMS as

an example. As MEMS functions over time, the growth of visible wear on rubbing surfaces between the gear and the pin joint often gives rise to the broken pin joint. The wear results from aging degradation and additional wear debris caused by some external shocks. Meanwhile, the external shocks can destroy the spring. Additionally, a micro-engine tends to deteriorate faster when its degradation level is greater than a certain threshold. For external shocks, their magnitudes are assumed to follow a normal distribution, where W ∼ N (µw , σw2 ). Combining preventive replacement, imperfect repair and corrective replacement, condition-based maintenance is also conducted for the MEMS based on the periodic inspection of the system state and the observation of the extent of visible wear. For a better exposition, the initially setting parameters are selected from the precedented researches with some modifications, and we add some parameters according to model assumptions shown in Table 2.

5.1

Reliability analysis

Since the complexity of Eq. (14), it is hard to calculate R(t) through the numerical integration, and we adopt Monte Carlo simulation to compute the system reliability. The algorithm is shown in the following pseudo-code:

19

Table 2. Parameters setting (a) Parameters setting for the system

Parameter H1 , H2 µc1 , σc1 µc2 , σc2 α1 , α2

(b) Parameters setting for random shocks

Value 1µm3 , 3µm3 2.5µm3 , 0.2µm3 3.5µm3 , 0.25µm3 0.06µm3 /Gpa, 0.065µm3 /Gpa

Parameter λ ∆ µw , σw D1 , D2 , D3

Value 50 0.1Gpa 1.2Gpa, 0.3Gpa 0.2Gpa, 1.2Gpa, 1.5Gpa

Algorithm 1 Monte Carlo simulation for R(t) Input: t Output: R(t) 1: set I = 0, M 2: for i ← 1 to M 3: do generate a realization {x(0), x(dt), · · · , x(T )} with the code provided in Appendix C 4: if x(t) < H2 5: then I = I + 1 6: end for 7: R(t) = I/M ;

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

R(t)

1

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

λ=0 λ=50 λ=100 λ=200

0 0

0.2

0.4

0.6

0.8

1

1.2

0

0.2

0.4

0.6

0.8

1

1.2

t

Fig. 9. Plot of the reliability function R(t)

Fig. 10. Sensitivity analysis of R(t) on λ

3

3

2.5

2.5

2

2

λ=50 λ=100 λ=200

Ratio

3.5

Ratio

3.5

1.5

1.5

1

1

0.5

0.5

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

t

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

Fig. 11. Plot of the ratio of the soft failure to the hard failure

Fig. 12. Sensitivity analysis of the ratio on λ

20

The system reliability function R(t) is shown in Fig. 9, and we can observe that the reliability function declines dramatically after t = 0.4, which suggests that the system transfers from the normal state to the weakened state, becoming more susceptible and having a larger risk to fail. Fig.10 presents the reliability curves with different occurrence rate λs, and it shows that R(t) shifts to the left with λ increasing. This result implies that the system would have a shorter lifetime if the external shock has a higher occurrence rate. Specifically, the green curve in Fig.10 is the reliability function of the system with no random shocks, and we can observe that random shocks shorten the system lifetime remarkably by almost 1/3 of the whole life cycle. For a better insight into the system failure process, Fig. 11 shows the ratio of the soft failure to the hard failure over time t. We can see that the curve increases drastically after t = 0.4 and then it flattens out as the ratio stops increasing. It suggests that system failures in the early life cycle are mainly caused by fatal shocks. However, with time going by, the soft failure accounts for about 75% of the system failures. In Fig. 12, as λ rises, the increasing of additional increments shortens the initial rise point of the curve. Meanwhile, random shocks with a high frequency decrease the maximum ratio of the soft failure to the hard failure because of more frequent fatal shocks.

5.2

The optimal CBM

Similarly, we compute CR(τ, H) with Monte Carlo simulation (see Algorithm 2) and the optimal values of decision variables for the proposed CBM strategy are determined through a simple numerical searching method in MATLAB. We set CI = 5, CR = 100, CP = 60, CIR = 20, τ1 = 0.3 and τ2 = 0.45. Note that the values of τ1 and τ2 depend on practical operation arrangements, and they could be altered in applications. Fig. 13 presents the surf and contour plots of CR with respect to H and τ . We can find that the minimum cost rate CR∗ is 7.464×10−3 at τ ∗ = 0.36 and H ∗ = 2.15. CR decreases before τ = τ ∗ while increases when τ > τ ∗ , which means too intensive inspections are unnecessary because of increasing inspection costs. In addition, at the optimal inspection period τ ∗ = 0.36, the reliability function begins to decrease drastically as shown in Fig. 9. Hence, it is reasonable to employ maintenance actions to keep a satisfactory

21

condition of the system, which can decrease the number of unexpected failures and improve system performances. Algorithm 2 Monte Carlo simulation for CR(τ, H) Input: τ, H Output: CR(τ, H) 1: set E(P ) = 0, E(C) = 0, CR , CP , CIR , CI 2: for i ← 1 to M 3: do generate a realization {x1 (0), x1 (dt), · · · , x1 (T )} with the code provided in Appendix D, where α = 0.06, D1 = 0.2, D2 = 1.2, D3 = 1.5, µc = 2.5, σc = 0.2. 4: set k1 = 0; k2 = 0 5: while x1 (k1 τ ) < H1 6: do k1 = k1 + 1 7: end 8: if x1 (k1 τ ) < H 9: then generate a realization {x2 (0), x2 (dt), · · · , x2 (T )} with the code pro10: vided in Appendix D, where α = 0.065, D1 = 0.1, D2 = 1.1 11: D3 = 1.4, µc = 3.5, σc = 0.25. 12: while x2 (k2 τ ) < H − H1 13: do k2 = k2 + 1; 14: end 15: if x2 (k2 τ ) < H2 − H1 16: then Cost = CP + CIR + (k1 + k2 )CI ; P = (k + k2 )τ 17: else Cost = CR + CIR + (k1 + k2 )CI ; P = (k1 + k2 )τ 18: elseif x1 (k1 τ ) < H2 19: then Cost = CP + k1 CI ; P = k1 τ 20: else Cost = CR + kCI ; P = k1 τ 21: E(P ) = P + E(P ), E(C) = Cost + E(C) 22: end for 23: CR(τ, H) = E(C)/E(P )

3

0.9

1.2

0.9

2.6

1

2.2 H

0.9

2

0.75

1.8

0.7

1 1.6 1.05

5

0.8

0.8 0.85

0.85

0.8

0.9 0.95

0.9 0.95

0.9

0.95

1

1.1

1.4 0.5

0.75

0.8

0.8

2.4

0.8

CR

5

0.8

1.1

0.85

0.8

2.8

1 1.5

1.2

2 0.35

2.5 0.3

τ

3

1 0.3

H

1.05

0.4

1.15

0.45

1 1.05 0.4

0.35 τ

(a) Surface plot of CR

(b) Contour plot of CR

Fig. 13. Surface and contour plots of CR with H and τ 22

1

1.1 1.15 0.45

5.3

Sensitivity analysis of CR∗

For realistic management operations, the maintenance cost could be cut down further through some adjustments in order to obtain a smaller CR∗ . So we conduct sensitivity analyses of CR∗ on CR , CP , and CIR shown in Table 3, Table 4, and Table 5 respectively. In the three tables, we list the 5%-20% reduction of maintenance costs including CR , CP and CIR , and the corresponding CR∗ improvements are computed. The results show that CR∗ is not very sensitive to either CR or CIR and every 5% change could only bring about 0.5%∼1% improvement of CR∗ . Comparatively, CR∗ is more sensitive to CP and every 5% cost reduction can improve around 3% of CR∗ , which indicates preventive maintenance accounting for a large proportion of the proposed CBM strategy. Therefore, it is wise to spend more efforts on reducing CP to obtain a better CR∗ . Note that these sensitivity analyses are subject to the initially setting parameters. Whereas in the managerial applications, the most appropriate decisions on reducing of CR , CP , and CIR depend on specific problems. Table 3. Sensitivity analysis of CR∗ on CR CR 100 95 90 80

Reduction of CR (%) 0 5 10 20

CR∗ (10−3 ) 7.464 7.421 7.374 7.222

Improvement of CR∗ (%) 0 0.576 1.206 3.242

Table 4. Sensitivity analysis of CR∗ on CP CP 60 57 54 48

Reduction of CP (%) 0 5 10 20

CR∗ (10−3 ) 7.464 7.248 7.005 6.582

Improvement of CR∗ (%) 0 2.894 6.149 11.817

Table 5. Sensitivity analysis of CR∗ on CIR CIR 20 19 18 16

Reduction of CIR (%) 0 5 10 20

CR∗ (10−3 ) 7.464 7.379 7.281 7.088

23

Improvement of CR∗ (%) 0 1.139 2.452 5.038

6

Conclusion and discussion In this work, a new reliability model with CBM for a binary-state deteriorating system subject

to external shocks is developed. Internal degradation of the system is modeled by a Wiener process, which has a larger shift and diffusion parameter when the system is in a weakened state. Random shocks are modeled by an HPP and classified into four zones, i.e., the SZ, the PDZ, the DZ, and the FZ. Among them, the FZ shock, which can devastate system structures directly, causes the hard failure. The PDZ and the DZ shocks contribute to the system deterioration and accelerate to the soft failure, and the SZ shock is considered harmless owing to the system strength. The dependency of zoned shock effects and the system state is discussed through an assumption of a reduction in classification thresholds of shock zones. To deal with failures of the binary-state deteriorating system, we propose a new CBM strategy with periodic inspections. Combining corrective replacement, preventive replacement, and imperfect repair, we optimize the average long-run total maintenance cost rate CR in the CBM strategy. Finally, a numerical example is presented to illustrate the applicability of the proposed reliability model with the optimal CBM strategy. While our main focus and contribution exist in the reliability and maintenance modeling parts, there are also some insights from the perspective of management. As shown in the numerical example, system failures in the early life cycle are mainly caused by fatal shocks. However, with time going by, the soft failure accounts for most system failures. Therefore, in the longrun operation, maintenance actions (e.g. lubrication, cleaning, etc.) which can decrease system wear is needed for industrial devices in order to keep a high efficiency of production. Apart from this, when the system is in the weakened state, the system reliability decreases at a higher speed. For this reason, it is suggested that maintenance actions should be implemented when the weakened state of operating systems is detected so as to reduce the unnecessary waste. For our proposed CBM, it has two characteristics: (1) the imperfect repair is included. (2) the maintenance action in CBM is implemented with a periodic schedule. Under this circumstance, the proposed CBM is easy to apply to industrial devices which need regular maintenance for the restrictions like geographical positions and the budget. A practical example is the maintenance offshore wind turbines. Even though their failures can be diagnosed continually, the immediate maintenance action may not be carried out if the operational condition is poor, e.g., extreme weather. Additionally, for some high-value products, repairing is a good way to decrease the maintenance cost and the proposed CBM combining with imperfect repair can reduce the long-

24

run maintenance cost rate to a great extent. Apart from what we have done in this work, there are some further issues that need to be addressed. In the numerical example part, we use Monte Carlo simulation to calculate the reliability function R(t) and the long run cost-rate CR(τ, H). Since we focus more on modeling rather than estimation, the uncertainty subject to data variability is not considered in this research. Nevertheless, it makes a great point for future studies. Another interesting topic is modeling multi-degradation processes with zoned shock effects. In most complex systems, there exist several degradation mechanisms like corrosion and fatigue. If degradation processes jointly interfere with environmental shocks, some different results may be obtained. In addition, we could add the δ-shock model to reconstruct zoned shock effects in both the magnitude dimension and the time dimension.

25

Appendix A As mentioned previously, the transition time Tw is the first passage time when the overall degradation level reaches H1 . It can be defined as: Tw = inf{t : X(t) ≥ H1 }. Since the overall degradation level is related to the shock process, the conditional CDF of Tw for a given number of the shocks in the DZ and PDZ is: F¯Tw |N2 (t),N3 (t) (t|n, m) = P(Tw > t| N2 (t) = n, N3 (t) = m) = P(X(t) < H1 | N2 (t) = n, N3 (t) = m) = P Xc1 (t) +

n X

(2)

Yj

j=0

+

m X

! (3)

Yj

< H1 |N2 (t) = n, N3 (t) = m

j=0





H1 − µc1 t − nµd1 − mµp1  . = Φ q 2 σc1 t + n2 σd2 + m2 σp2 Based on Eqs. (2) and (16), we have the CDF of Tw as:

FTw (t) = 1 −

∞ X ∞ X





H1 − µc1 t − nµd1 − mµp1  e−λp2 t (λp2 t)n e−λp3 t (λp3 t)m Φ q , n! m! 2 σc1 t + n2 σd2 + m2 σp2 n=0 m=0

where Φ(·) is the CDF of the standard normal distribution.

Appendix B Similar to FTw (·), we derive the CDF FTH (·) of TH . As mentioned previously, TH is the first passage time when the overall degradation level reaches H with a given initial degradation level H1 at time jτ . It can be defined as: TH = jτ + inf {t : X(jτ + t) − X(jτ ) ≥ H − H1 } .

26

Denote TH0 = TH − jτ , and for a given number of N200 (t) and N300 (t), the CDF of TH0 is: F¯TH0 |N200 (t),N300 (t) (t|n, m) = P(TH0 > t|N200 (t) = n, N300 (t) = m) = P ∆Xc2 (t) +

n X

(2)

Yj

j=0

 = Φ

+

m X

! (3)

Yj

< H − H1

j=0



H − H1 − µc2 t − nµd2 − mµp2  q , 2 σc2 t + n2 σd2 + m2 σp2

where ∆Xc2 (t) is the change of continuous degradation from time jτ to jτ + t. N200 (t) and N300 (t) are the counting numbers of random shocks in the PDZ and the DZ from time jτ to time jτ + t for a system in the weakened state respectively. Then, the unconditional CDF of TH is derived as: FTH (t) = 1 −

∞ X ∞ X

F¯TH0 |N200 (t+jτ ),N300 (t+jτ ) (t + jτ |n, m) · P{N200 (t + jτ ) = n} · P{N300 (t + jτ ) = m}

n=0 m=0 ∞ X ∞ X

0

0

e−λp2 (t+jτ ) [λp02 (t + jτ )]n e−λp3 (t+jτ ) [λp03 (t + jτ )]m =1− n! m! n=0 m=0   H − H1 − µc2 (t + jτ ) − nµd2 − mµp2  q Φ . 2 2 2 2 2 σc2 (t + jτ ) + n σd + m σp

Appendix C The MATLAB code:

1

function [Y]=Algorithm1(lambda,dt,T)

2

alpha 1=0.06;alpha 2=0.065;H 1=1;H 2=3;mu w=1.2;sigma2 w=0.09;

3

D 1=0.2;D 2=1.2;D 3=1.5;mu c1=2.5;sigma c1=0.2;

4

mu c2=3.5,sigma c2=0.25;D 12=0.1;D 22=1.1;D 32=1.4;

5

A=[];

6

l=exprnd(1/lambda);

7

while l
8

A=[A,l];

9

l=l+exprnd(1/lambda);

10

end

11

Y=[];k=0;s=0;a=0;c=0;

12

while s<1;

13

if any(fix(A*10000)/10000==k*dt)==1;%see if the shock arrives at the ...

27

system. sign=1;

14

else

15

sign=0;

16 17

end

18

w=normrnd(mu w,sigma2 w);%add compound HPP to Winere process.

19

if w>D 3; d=H 2 *sign+c;

20

else if w>D 2;

21

d=alpha 1*w*sign+c;

22

else if w>D 1

23

d=alpha 1 *w*(w-D 1)/(D 2-D 1)*sign+c;

24

end

25

end

26 27

end

28

b=a+sqrt(dt)*randn;

29

s=mu c1 *(k*dt)+sigma c1*b+d;

30

Y=[Y,s];

31

k=k+1;

32

a=b;

33

c=d;

34

end

35

f=0;g=0;h=0;

36

for i=k*dt:dt:T;

37

if any(fix(A*10000)/10000==i)==1;%see if the shock arrivels at the system. sign=1;

38 39

else sign=0;

40 41

end

42

w=normrnd(mu w,sigma2 w);%add compound HPP to Winere process.

43

if w>D 32; h=H 2 *sign+g;

44 45

else if w>D 22; h=alpha 1*w*sign+g;

46

else if w>D 12

47

h=alpha 2*w*(w-D 12)/(D 22-D 12)*sign+g;

48

end

49

end

50 51

end

52

e=f+sqrt(dt)*randn;

28

53

m=mu c2 *(i-k*dt)+sigma c2*e+H 1+h;

54

Y=[Y,m];

55

f=e;

56

g=h;

57

end

58

end

Appendix D The MATLAB code:

1

function [Y]=Algorithm2(alpha,D 1,D 2,D 3,mu c,sigma c,dt,T)

2

mu w=1.2;sigma2 w=0.09;lambda=50;H 2=3;

3

A=[];

4

l=exprnd(1/lambda);

5

while(l
6

A=[A,l];

7

l=l+exprnd(1/lambda);

8

end

9

Y=[];k=0;s=0;a=0;c=0;

10 11

for t=0:dt:T; if any(fix(A/dt)==k)==1;%see if the shock arrives at the system. sign=1;

12 13

else sign=0;

14 15

end

16

w=normrnd(mu w,sigma2 w);%add compound HPP to Winere process.

17

if w>D 3; d=H 2 *sign+c;

18 19

else if w>D 2; d=alpha*w*sign+c;

20

else if w>D 1

21

d=alpha*w*(w-D 1)/(D 2-D 1)*sign+c;

22

end

23

end

24 25

end

26

b=a+sqrt(dt)*randn;

27

s=mu c *(k*dt)+sigma c *b+d;

29

28

Y=[Y,s];

29

k=k+1;

30

a=b;

31

c=d;

32

end

33

end

References Alaswad, S. and Xiang, Y. (2017). A review on condition-based maintenance optimization models for stochastically deteriorating system. Reliability Engineering & System Safety, 157:54–63. Bae, S. J. and Kvam, P. H. (2006). A change-point analysis for modeling incomplete burn-in for light displays. IIE Transactions, 38(6):489–498. Castro, I. T., Caball, N. C., and Prez, C. J. (2015). A condition-based maintenance for a system subject to multiple degradation processes and external shocks. International Journal of Systems Science, 46(9):1692–1704. Do, P. and B´erenguer, C. (2012). Condition-based maintenance with imperfect preventive repairs for a deteriorating production system. Quality & Reliability Engineering International, 28(6):624–633. Do, P., Voisin, A., Levrat, E., and Iung, B. (2015). A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions. Reliability Engineering & System Safety, 133:22–32. Finkelstein, M. S. and Zarudnij, V. I. (2001). A shock process with a non-cumulative damage. Reliability Engineering & System Safety, 71(1):103–107. Huang, Z., Xu, Z., Wang, W., and Sun, Y. (2015). Remaining useful life prediction for a nonlinear heterogeneous wiener process model with an adaptive drift. IEEE Transactions on Reliability, 64(2):687–700. Huynh, K. T., Castro, I. T., Barros, A., and B´erenguer, C. (2012). Modeling age-based maintenance strategies with minimal repairs for systems subject to competing failure modes due to degradation and shocks. European Journal of Operational Research, 218(1):140–151.

30

Jiang, L., Feng, Q., and Coit, D. W. (2015). Modeling zoned shock effects on stochastic degradation in dependent failure processes. IIE Transactions, 47(5):460–470. Keizer, M. C. O., Flapper, S. D. P., and Teunter, R. H. (2017). Condition-based maintenance policies for systems with multiple dependent components: A review. European Journal of Operational Research, 261(2):405–420. Lawless, J. and Crowder, M. (2004). Covariates and random effects in a gamma process model with application to degradation and failure. Lifetime Data Analysis, 10(3):213–227. Lehmann, A. (2009). Joint modeling of degradation and failure time data. Journal of Statistical Planning and Inference, 139(5):1693–1706. Liu, B., Xie, M., Xu, Z., and Kuo, W. (2016). An imperfect maintenance policy for missionoriented systems subject to degradation and external shocks. Computers & Industrial Engineering, 102:21–32. Nakagawa, T. (2007). Shock and damage models in reliability theory. Springer Science & Business Media. Peng, C.-Y. (2015). Inverse gaussian processes with random effects and explanatory variables for degradation data. Technometrics, 57(1):100–111. Peng, H., Feng, Q., and Coit, D. W. (2010). Reliability and maintenance modeling for systems subject to multiple dependent competing failure processes. IIE transactions, 43(1):12–22. Rafiee, K., Feng, Q., and Coit, D. W. (2014). Reliability modeling for dependent competing failure processes with changing degradation rate. IIE transactions, 46(5):483–496. Rafiee, K., Feng, Q., and Coit, D. W. (2015). Condition-based maintenance for repairable deteriorating systems subject to a generalized mixed shock model. IEEE Transactions on Reliability, 64(4):1164–1174. Ross, S. M. (2014). Introduction to probability models. Academic Press. Wang, X. and Xu, D. (2010). An inverse gaussian process model for degradation data. Technometrics, 52(2):188–197.

31

Yang, L., Ma, X., Peng, R., Zhai, Q., and Zhao, Y. (2017a). A preventive maintenance policy based on dependent two-stage deterioration and external shocks. Reliability Engineering & System Safety, 160:201–211. Yang, L., Ma, X., and Zhao, Y. (2017b). A condition-based maintenance model for a three-state system subject to degradation and environmental shocks. Computers & Industrial Engineering, 105:210–222. Ye, Z.-S., Xie, M., Tang, L.-C., and Shen, Y. (2012). Degradation-based burn-in planning under competing risks. Technometrics, 54(2):159–168. Zhang, J., Huang, X., Fang, Y., Zhou, J., Zhang, H., and Li, J. (2016). Optimal inspection-based preventive maintenance policy for three-state mechanical components under competing failure modes. Reliability Engineering & System Safety, 152:95–103. Zhang, Z., Hu, C., Si, X., Zhang, J., and Shi, Q. (2017). A prognostic approach for systems subject to wiener degradation process with cumulative-type random shocks. In Data Driven Control and Learning Systems (DDCLS), 2017 6th, pages 694–698. IEEE. Zhang, Z., Si, X., Hu, C., and Lei, Y. (2018). Degradation data analysis and remaining useful life estimation: A review on wiener-process-based methods. European Journal of Operational Research.

32

Highlights (1) This study considers a deteriorating system with a two-phase degradation path based on the overall degradation level. (2) Shock effects are comprehensively incorporated in the deteriorating system, and they are assumed to change with the system state. (3) A new condition-based maintenance strategy is proposed integrating preventive replacement, corrective replacement and imperfect repair. (4) The average long-run maintenance cost rate is derived to optimize the condition-based maintenance strategy.

1

Acknowledgement This research was supported by National Natural Science Foundation of China (Grants No. 71472132, 71532008, &71802145).

1