Reliability Evaluation and Optimization of a System with Mixed Run Shock

Abstract

In this paper, we investigate a wear and mixed shock model in which the system can fail due to internal aging or external shocks. The lifetime of the system, due to internal wear, follows continuous phase-type (PH) distributions. The external random shocks arrive at the system according to a PH renewal process. The system will fail when the internal failure occurs or $k_1$ consecutive external shocks, the size of at least $d_1$ or $k_2$ consecutive external shocks the size of at least $d_2$ occur, where $d_1<d_2$, $k_1>k_2$. The failed system can be repaired immediately, and the repair times of the system are governed by continuous PH distributions. The system can be replaced by a new and identical one based on a bivariate replacement policy $(L,N)$. The long-run average profit rate for the system is obtained by employing the closure property of the PH distribution. Finally, a numerical example is also given to determine the optimal replacement policy.

1. Introduction

In engineering practice, systems usually operate in complex environments [1,2]. In addition to internal wear degradation, systems are frequently impacted by the external environment, such as voltage fluctuations, temperature and humidity change, etc. [3,4]. The external factors can be regarded as shocks attacking the systems [5]. For example, a system may fail due to the change in its internal structure or the fluctuation of the external environmental conditions, such as strength, temperature, voltage, and vibrations. The existing shock models are generally classified into four categories, namely, extreme, cumulative, run, and $\delta$-shock models. The extreme shock model is one of the most widely used shock models, in which the system will fail when the size of any shock is beyond a specified threshold value. Gut and Hüsler [6] derived moment relations and asymptotic distributions of the time to failure under the extreme shock model. Eryilmaz and Kan [7] investigated reliability and optimal replacement policy for an extreme shock model with a change point. In the run shock model, failure occurs when there is a run of $n$ consecutive shocks that are greater than a threshold level. In recent years, researchers also proposed some new mixed shock models. Lorvand et al. [8] studied the life behavior of a shock model, in which the system fails when k out of interarrival times between two successive shocks is less than δ, or the magnitude of the shock is larger than the other critical level, say $\gamma$. Parvardeh and Balakrishnan [9] studied the system lifetime mixed $\delta$-shock models, in which the system can fail when the time between two successive shocks is less than a critical threshold $\delta$, or the magnitude of the shock (the cumulative magnitude of shocks) is larger than another critical threshold $\gamma$. Mallor and Omey [10] proposed systems that fail when k consecutive shocks with critical magnitude (e.g., above or below a certain critical level) occur. Eryilmaz and Tekin [11] evaluated system reliability under a mixed shock model, in which the system under concern failed upon the occurrence of $k$ consecutive shocks the size of at least $d_1$ or a single large shock the size of at least $d_2$.

In addition to external random shocks, internal wear degradation is another important cause of system failure. Moreover, the failure caused by external shocks and internal wear degradation are often competitive with each other. Reliability modeling for systems with multiple dependent competing failure processes (MDCFP) attracted more attention in recent years. Peng et al. [12] proposed reliability and maintenance models for systems subject to MDCFP with a changing, dependent failure threshold. Rafiee et al. [13] presented reliability models for devices subject to dependent competing failure processes of degradation and random shocks with a changing degradation rate, in which four different shock patterns that can increase the degradation rate were considered. Jiang et al. [14] studied a system that experiences two dependent competing failure processes. Qiu and Cui [15] evaluated system reliability performance based on a dependent two-stage failure process with competing failures. Liu et al. [16] proposed a novel system subject to MDCFP under chance theory. In the previous study of MDCFP, the lifetime of components and the interarrival time of shocks were usually assumed to be exponentially distributed. The exponential distribution models the behavior of components that fail at a constant rate, regardless of the accumulated age. Although this property greatly simplifies analysis, it makes the distribution inappropriate for some real-world scenarios. To overcome this issue, PH distributions were introduced by Marcel Netus in 1975 (Neuts [17]) and are used widely in science and engineering. The set of PH-distributions has a number of useful properties: (i) They can approximate any probability distribution with a nonnegative support, which leads to high applicability. (ii) They are associated with finite state Markov chains and can be utilized to introduce Markov processes for stochastic systems. (iii) Finally, the set of distributions is closed under a number of basic operations, including addition, min, and max, which leads to exact and closed form solutions. In recent years, the PH distribution was also used in the reliability field by several scholars. Kim and Kim [18] presented reliability models for a nonrepairable system with heterogeneous components having phase-type time-to-failure distributions. Liu et al. [19] developed a cold standby repairable system by adopting working vacations, vacation interruption policies, and PH distributions. Wen et al. [20] investigated a multiple warm standby delta shock system by considering the multiple vacations policy and PH distributions. Montoro-Cazorla and Pérez-Ocón [21,22,23] proposed different shock and wear models by using PH distributions and Markovian arrival processes (MAPs). Considering that the number of shocks the device can support is limited, Montoro-Cazorla et al. [24] presented a device submitted to shocks and wear by using PH distributions. Yu et al. [25] studied the optimal order replacement policy for a repairable system by employing PH distributions, in which a bivariate maintenance policy $(K,N)$ was adopted. Considering the perspective of rational use of human resources, Yu et al. [26] investigated the reliability for a phase-type geometric process repair model with spare device procurement and a repairman’s multiple vacations. Eryilmaz [27] computed the optimal replacement time and mean residual lifetime of a system under particular class of reliability shock models by employing the closure property of the PH distribution.

In this paper, a multiple competing failure geometric process repair system model is proposed by combing PH distributions and geometric processes. The lifetime of the system, due to internal wear, the repair time for repair of the failure system, replacement time for replacement of the failure system, and the interarrival time of the shocks are governed by different PH distributions, respectively. In addition, the consecutive operating time after failure is getting shorter and shorter, and the consecutive maintenance time after failure will become longer and longer by employing the well-known geometric processes. In addition to internal wear degradation, external random shocks are another cause of system failure. The external shocks are a new mixed shock model in which systems will fail when $k_1$ consecutive external shocks the size of at least $d_1$ or $k_2$, consecutive external shocks the size of at least $d_2$ occurs, where $d_1<d_2$, $k_1>k_2$. Clearly, if $d_1\to \infty$ or $d_2\to \infty$, then the new mixed shock model is the same with the run shock model. If $k_1=1,d_2\to \infty$ or $k_2=1,d_1\to \infty$, then the new mixed shock is the same with the extreme shock model, i.e., the extreme and run shock model are special cases of the new mixed shock model. Under this new shock model, a small number of consecutive shocks of large damage sizes, or a large number of consecutive shocks of small damage sizes can cause system failure. Moreover, the internal wear and external random shocks are competing with each other. Furthermore, a bivariate replacement policy $(L,N)$ is adopted. Here, the symbol $L$ means that the system is replaced when the total working time of the system exceeds $L$ time units, while $N$ means that the system is replaced at the time of the $N\mathrm{th}$ failure, whichever occurs first. Compared with previous studies, the innovations of this paper are summarized as follows: (1) a multiple competing failure geometric process repair system model is proposed; (2) a new mixed run shock pattern is introduced into the model; (3) the bivariate replacement policy $(L,N)$ is adopted; and (4) the explicit express of the long-run average profit rate function under replacement policy $(L,N)$ is obtained.

In fact, the mixed run shock pattern was put forward for future work by Eryilmaz and Tekin [11]. However, the problem was not solved by researchers until now. The mixed run shock pattern is considered for the multiple competing failure geometric process repair system model in this paper.

The rest of this paper is organized as follows: In Section 2, the model assumption is described, and some related preliminaries are also presented. The mean time between two consecutive replacements is discussed in Section 3. In Section 4, the explicit express of the long-run average profit rate function under replacement policy $(L,N)$ is obtained. A numerical example to illustrate the proposed model is given in Section 5. Finally, the conclusions are summarized in Section 6.

2. The Model Assumption and Preliminaries

The detailed assumptions of the system with mixed run shock are given as follows:

Assumption 1.

Consider a system that is subject to internal wear degradation and a sequence of external shocks over time. At time$t=0$, a new system is put into operation. It is obvious that the system can fail in two ways: the internal wear degradation due to the materiel fatigue and ageing and the external failure caused by the frequency of shocks.

Assumption 2.

A bivariate replacement policy$(L,N)$based on the number of failures and total working time of the system, is adopted, under which the system is replaced at the time of the$N\mathrm{th}$failure or the total working time of the system exceeds$L$time units, whichever occurs first.

Assumption 3.

The time interval between starting the operation of a new system and the completion of the first repair or replacement is defined as the first sub-cycle of the system. The time interval between completion of the $(n-1)\mathrm{th}$repair and completion of the$n\mathrm{th}$repair or the first replacement is defined as the$n\mathrm{th}$sub-cycle of the system, $n=1,2,\dots ,N-1$. Additionally, the time interval from the end of the$(N-1)\mathrm{th}$repair to the completion of the replacement is called the $N\mathrm{th}$sub-cycle of the system. The time internal between two consecutive replacements is called a full cycle. Furthermore, we assume that the random replacement time$R$for replacement of the failure system has a continuous PH distribution with representation$P{H}_c(\delta ,\mathit{H})$of order$m_1$.

Assumption 4.

Let$X_i^n\hspace{0.17em}(n=1,2,\dots ,N,i=1,2,\dots )$be the interarrival time between the$(i-1)\mathrm{th}$and the$i\mathrm{th}$external shocks in the$n\mathrm{th}$sub-cycle. We assume that the interarrival times$X_i^n$are independent identically distributed random variables, and follow a common continuous PH distribution, represented as$P{H}_c(\beta ,S)$with order$s$. Define$Y_i^n(n=1,2,\dots ,N,i=1,2,\dots )$to be the magnitude of the$i\mathrm{th}$shock in the$n\mathrm{th}$sub-cycle, and suppose$Y_i^n$are independent identically distributed random variables. Let us fix the two threshold values,$d_1$and$d_2$, such that$d_1<d_2$, and the system fails due to external shocks when$k_1$consecutive shocks the size of at least$d_1$or$k_2$consecutive shocks the size of at least$d_2$occurs, where$k_1>k_2$.

Assumption 5.

Let$W_n(n=1,2,\dots ,N)$be the lifetime of the system due to internal wear degradation in the$n\mathrm{th}$sub-cycle. The random variable$W_n$follows a continuous PH distribution with representation$P{H}_c(\alpha ,a^{n-1}T)$of order$m_2$, where$a$is a real constant and$a\ge 1$. Obviously, the sequence$\{W_n,n=1,2,\dots ,N\}$is a decreasing geometric process with ratio$a$.

Assumption 6.

Let$Z_n(n=1,2,\dots ,N-1)$be the repair time after the$n\mathrm{th}$failure, and$Z_n$follows a continuous PH distribution with representation$P{H}_c(\gamma ,b^{n-1}G)$of order$m_3$, where$b$is a real constant and$0<b<1$. Then the sequence$\{Y_n,n=1,2,\dots ,N-1\}$forms an increasing geometric process with ratio$b$.

Assumption 7.

The random variables$R$, $X_i^n$, $Y_i^n$, $W_n$, and$Z_n$are independent of each other.

Let $N_n$ be the number of shocks until $k_1$ consecutive shocks of a size above or equal to $d_1$ or $k_2$ consecutive shocks the size of at least $d_2$ in the $n\mathrm{th}$ sub-cycle, then the lifetime of the system due to external shocks in the $n\mathrm{th}$ sub-cycle can be denoted as

$$T_n={\displaystyle \sum _{i=1}^{N_n}{X}_i^n}$$

(1)

For better understanding of the proposed model, Figure 1 presents a possible realization of external shocks of the system. According to the assumption of the proposed model, if $k_1=3,k_2=2$, the system will fail when three consecutive shocks the size of at least $d_1$ or two consecutive shocks the size of at least $d_2$ occurs. From Figure 1, it is obvious that the system fails after the fifth shock because of the magnitudes of the third, fourth, and fifth shocks, which are greater than $d_1$. Hence the system’s lifetime is $T_n=X_1^n+\cdots +X_5^n$. If $k_1=5,k_2=2$, then the system failure occurs due to the magnitudes of the eighth and ninth shocks, which are greater than $d_2$. In this case, the system’s lifetime is $T_n=X_1^n+\cdots +X_9^n$.

It is obvious that the system fails when the internal wear degradation process enters the absorb state or $k_1$ consecutive shocks the size of at least $d_1$ or $k_2$ consecutive shocks the size of at least $d_2$ occurs, whichever occurs first. Therefore, the lifetime of the system due to internal wear degradation and external shocks in the $n\mathrm{th}$ sub-cycle can be represented as

$$O_n=\min (W_n,T_n)$$

(2)

Theorem 1.

Let$p_1=P\{d_1\le Y_i^n<d_2\}$and$p_2=P\{Y_i^n\ge d_2\}$for$i=1,2,\dots$, then the random variable$N_n$follows a discrete PH-distribution represented as$P{H}_d(\epsilon ,\mathit{Q})$with order$k_1+k_2-1$, where$\epsilon ={(1,0,\dots ,0)}_{{}_{1\times [(k_1+k_2-1)\times (k_1+k_2-1)]}}$,

$$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{ccccccccc}0& \hspace{1em}\hspace{0.17em}1& \hspace{0.17em}11& \cdots & \underset{k_1-1}{\underbrace{11\cdots 1}}& 2& 22& \cdots & \underset{k_2-1}{\underbrace{22\cdots 2}}\end{array}\\ \mathit{Q}=\begin{array}{c}\\ 0\\ 1\\ 11\\ ⋮\\ \underset{k_1-1}{\underbrace{11\cdots 1}}\\ 2\\ 22\\ ⋮\\ \underset{k_2-1}{\underbrace{22\cdots 2}}\end{array}{\left(\begin{array}{ccccccccc}1-p_1-p_2& p_1& 0& \cdots & \hspace{0.17em}0& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_2& 0& \cdots & \hspace{1em}0\\ 1-p_1-p_2& 0& p_1& \cdots & \hspace{0.17em}0& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_2& 0& \cdots & \hspace{1em}0\\ 1-p_1-p_2& 0& 0& \cdots & \hspace{0.17em}0& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_2& 0& \cdots & \hspace{1em}0\\ ⋮& ⋮& ⋮& \ddots & ⋮& \hspace{1em}⋮& ⋮& \cdots & \hspace{1em}0\\ 1-p_1-p_2& 0& 0& \cdots & 0& \hspace{1em}0& 0& \cdots & \hspace{1em}0\\ 1-p_1-p_2& p_1& 0& \cdots & 0& \hspace{1em}0& p_2& \cdots & \hspace{1em}0\\ 1-p_1-p_2& p_1& 0& \cdots & 0& \hspace{1em}0& 0& \cdots & \hspace{1em}0\\ ⋮& ⋮& ⋮& \ddots & ⋮& \hspace{1em}⋮& ⋮& \ddots & \hspace{1em}⋮\\ 1-p_1-p_2& p_1& 0& \cdots & 0& \hspace{1em}0& 0& \cdots & \hspace{1em}0\end{array}\right)}_{(k_1+k_2-1)\times (k_1+k_2-1)}\end{array}$$

(3)

and$\mathit{u}={(0,0,\dots ,0,p_1+p_2,0,\cdots ,0,p_2)}_{1\times (k_1+k_2-1)}$.

Proof. 

Let us define a sequence of trials $I_i,i=1,2,\dots$ such that

$$I_i=\{\begin{array}{l}0,\hspace{1em}\mathrm{if}\hspace{0.17em}{Y}_i^n<d_1,\\ 1,\hspace{1em}\mathrm{if}\hspace{0.17em}{d}_1\le Y_i^n<d_2,\\ 2,\hspace{1em}\mathrm{if}\hspace{0.17em}{Y}_i^n\ge d_2.\end{array}$$

(4)

and the random variable $N_n$ is the time to absorption for a Markov chain defined by the sequence $I_i,i=1,2,\dots$. Then, the corresponding Markov chain has $k_1+k_2-1$ transient states, which can be denoted as $\{0,1,11,\dots ,\underset{k_1-1}{\underbrace{11\cdots 1}},2,22,\dots ,\underset{k_2-1}{\underbrace{22\cdots 2}}\}$, and three absorbing states $\{\underset{k_1}{\underbrace{11\cdots 1}},\underset{k_1-1}{\underbrace{11\cdots 1}}2,\underset{k_2}{\underbrace{22\cdots 2}}\}$. The transition probability matrix of the Markov chain between transient states is given by $Q$, and the transition probability vector from transient states to absorbing states is given by $u$. The proof is completed. □

Notice: According to the definition of discrete PH distributions (He [28]), a discrete PH distribution is the distribution of the time to absorbing in an absorbing Markov chain, and its probability mass function is as follows:

$$P\{N=n\}=\epsilon Q^{n-1}u,\hspace{0.17em}n=1,2,\dots .$$

Since the $k_1$-th element in the first row of matrix $Q^{n-1}$ denotes the probability that in consecutive $n-1$ shocks, the sizes of the last $k_1-1$ consecutive shocks are at least $d_1$ but not more than $d_2$, or the sizes of the last $k_2-1$ consecutive shocks are at least $d_2$, and the $k_1+k_2-1$-th element in the first row of matrix $Q^{n-1}$ denotes the probability that in consecutive $n-1$ shocks, the sizes of the last $k_2-1$ consecutive shocks are at least $d_2$, then $\epsilon Q^{n-1}u$ denotes the probability that in consecutive $n$ shocks, the sizes of the last $k_1$ consecutive shocks are at least $d_1$ but not more than $d_2$, or the sizes of the last $k_2$ consecutive shocks are at least $d_2$, that is, the Markov chain enters absorbing states. Similar methodology to proof Theorem 1 can be found in Lemma 1 of reference [11].

Theorem 2.

If the interarrival times$X_i^n,i=1,2,\dots$are independent identically distributed random variables, and$X_i^n\sim P{H}_c(\beta ,S)$with order$s$, and independently the number of shocks until system failure due to external shocks$N_n\sim P{H}_d(\epsilon ,Q)$with order$k_1+k_2-1$. Then, the system’s lifetime, due to external shocks in the$n\mathrm{th}$sub-cycle, follows a continuous PH distribution of the order$s(k_1+k_2-1)$with representation

$$T_n={\displaystyle \sum _{i=1}^{N_n}{X}_i^n}\sim P{H}_c(\beta \otimes \epsilon ,S\otimes I+(S^0\beta )\otimes Q)$$

(5)

Proof. 

According to the well-known closure property of PH distributions (He [28]), the distribution of $T_n$ can be immediately obtained. The proof is completed. □

Theorem 3.

If the lifetime of the system, due to internal wear degradation, is$W_n\sim P{H}_c(\alpha ,a^{n-1}T)$with order$m_2$, and independently the lifetime of the system due to external shocks is$T_n\sim P{H}_c(\beta \otimes \epsilon ,S\otimes I+(S^0\beta )\otimes Q)$with order$s(k_1+k_2-1)$. Then, the system’s lifetime$O_n$, due to internal wear degradation or external shocks in the$n\mathrm{th}$sub-cycle, follows a continuous PH distribution of order$m_{2}s(k_1+k_2-1)$with representation

$$O_n\sim P{H}_c(\alpha \otimes (\beta \otimes \epsilon ),a^{n-1}T\oplus (S\otimes I+(S^0\beta )\otimes Q))$$

(6)

Proof. 

The distribution of the system’s lifetime $O_n$ due to internal wear degradation or external shocks in the $n\mathrm{th}$ sub-cycle can be can be immediately obtained by Proposition 1.3.1 of He [28]. The proof is completed. □

According to the closure properties of the PH distributions, the following results can be obtained:

(a)

The random variable ${\displaystyle \sum _{i=1}^n{O}_i},(n=2,3,\dots ,N)$ follows a continuous PH distribution of order $n{m}_{2}s(k_1+k_2-1)$ with representation $P{H}_c({\phi }_n,{\Delta }_n)$, where

$${\phi }_n=(\alpha \otimes (\beta \otimes \epsilon ),{\mathit{O}}_{1\times (n-1)m_{2}s(k_1+k_2-1)})$$

(7)

$${\Delta }_n=\left(\begin{array}{ccccc}{\Delta }_{11}& {\Delta }_{12}& & & \\ & {\Delta }_{22}& {\Delta }_{23}& & \\ & & \ddots & \ddots & \\ & & & {\Delta }_{(n-1)\times (n-1)}& {\Delta }_{(n-1)\times n}\\ & & & & {\Delta }_{nn}\end{array}\right),$$

(8)

$${\Delta }_{kk}=a^{k-1}T\oplus (S\otimes I+(S^0\beta )\otimes Q),k=1,2,\dots ,n,$$

(9)

$${\Delta }_{k\times (k+1)}=-(a^{k-1}T\oplus (S\otimes I+(S^0\beta )\otimes Q))e(\alpha \otimes (\beta \otimes \epsilon )),k=1,2,\dots ,n-1.$$

(10)

(b)

The random variable ${\displaystyle \sum _{n=1}^N{Z}_n}$ follows a continuous PH distribution of order $N{m}_3$ with representation $P{H}_c(\tilde{\phi },\tilde{\Delta })$, where $\tilde{\phi }=(\gamma ,{\mathit{O}}_{(N-1)m_3})$,

$$\tilde{\Delta }=\left(\begin{array}{ccccc}G& -Ge\gamma & & & \\ & bG& -bGe\gamma & & \\ & & \ddots & \ddots & \\ & & & b^{N-2}G& -b^{N-2}Ge\gamma \\ & & & & b^{N-1}G\end{array}\right)$$

(11)

3. Mean Time between Replacement (MTBR)

The time length of two consecutive replacements is equal to the sum of working time, repair time, and replacement time in a full cycle, thus the MTBR under $(L,N)$ policy can be represented as follows

$${\mathrm{MTBR}}_{(L,N)}=E[\min \{{\displaystyle \sum _{n=1}^N{O}_n},L\}]+E[{({\displaystyle \sum _{n=1}^N{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}}+{({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}>L\}}]+E[R],$$

(12)

where ${\chi }_A$ is an indicator function such that

$${\chi }_A=\{\begin{array}{l}1,\hspace{0.17em}\mathrm{if}\hspace{0.17em}\mathrm{event}\hspace{0.17em}A\hspace{0.17em}\mathrm{occurs},\\ 0,\hspace{0.17em}\mathrm{if}\hspace{0.17em}\mathrm{event}\hspace{0.17em}\mathrm{doesnot}\hspace{0.17em}A\hspace{0.17em}\mathrm{occur}.\end{array}$$

and $M(L)$ denotes the number of completed sub-cycles by time $L$.

Remark 1.

When$M(L)=0$, the system operates in time interval$[0,L]$, then is replaced at time$t=L$by a new and identical one; it is a very simple case of our proposed model, thus it is not considered in our present model.

In the following, we compute the mean working time, mean repair time, and mean replacement time in a full cycle.

The mean working time in a full cycle is

$$\begin{array}{l}E[\min \{{\displaystyle \sum _{n=1}^N{O}_n},L\}]={\displaystyle {\int }_0^{L}x{\phi }_N\exp ({\Delta }_{N}x){\Delta }^{0}dx}+{\displaystyle {\int }_L^{+\infty }L{\phi }_N\exp ({\Delta }_{N}x){\Delta }_N^{0}dx}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum _{n=0}^{\infty }\frac{(n+1)}{(n+2)!}}{\phi }_N{\Delta }_N^{-2}{({\Delta }_{N}L)}^{n+2}{\Delta }_N^0+L{\phi }_N\exp ({\Delta }_{N}L)e.\end{array}$$

(13)

The mean repair time in a full cycle can be denoted as

$$E[{({\displaystyle \sum _{n=1}^N{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}}+{({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}>L\}}],$$

(14)

The first term summation in Equation (14) can be computed by

$$\begin{array}{l}E[{({\displaystyle \sum _{n=1}^N{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}}]=E[{\displaystyle \sum _{n=1}^N{Z}_n}|{\displaystyle \sum _{n=1}^N{O}_n}\le L]\cdot P\{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}={\displaystyle \sum _{n=1}^N(-\gamma b^{n-1}{G}^{-1}e)\cdot }(1-{\phi }_N\exp ({\Delta }_{N}L)e).\end{array}$$

(15)

The second term summation in Equation (14) is given by

$$\begin{array}{l}E[{({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}>L\}}]=E[({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})|{\displaystyle \sum _{n=1}^N{O}_n}>L]\cdot P\{{\displaystyle \sum _{n=1}^N{O}_n}>L\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=E[E[({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})|M(L)]]\cdot {\phi }_N\exp ({\Delta }_{N}L)e\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=E[{\displaystyle \sum _{n=1}^{N-1}{Z}_n}|M(L)=n]\cdot P\{M(L)=n\}\cdot {\phi }_N\exp ({\Delta }_{N}L)e\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum _{n=1}^{N-1}(-\gamma b^{n-1}{G}^{-1}e)}\cdot P\{M(L)=n\}\cdot {\phi }_N\exp ({\Delta }_{N}L)e,\end{array}$$

(16)

On the other hand, we have

$$P\{M(L)=n\}=P\{{\displaystyle \sum _{k=1}^n{O}_n}<L\le {\displaystyle \sum _{k=1}^{n+1}{O}_n}\}={\phi }_{n+1}\exp ({\Delta }_{n+1}L)e-{\phi }_n\exp ({\Delta }_{n}L)e$$

(17)

Based on Equations (15)–(17), the mean repair time in a full cycle can be obtained

$$\begin{array}{l}E[{({\displaystyle \sum _{n=1}^N{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}}+{({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}>L\}}]={\displaystyle \sum _{n=1}^N(-\gamma b^{n-1}{G}^{-1}e)\cdot }(1-{\phi }_N\exp ({\Delta }_{N}L)e)\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^{N-1}(-\gamma b^{n-1}{G}^{-1}e)}\cdot ({\phi }_{n+1}\exp ({\Delta }_{n+1}L)e-{\phi }_n\exp ({\Delta }_{n}L)e)\cdot {\phi }_N\exp ({\Delta }_{N}L)e.\end{array}$$

(18)

From the model assumptions we obtain

$$E[R]=-\delta {\mathit{H}}^{-1}{e}_{m_1}$$

(19)

Substituting Equations (13), (18) and (19) into (14), we obtain

$$\begin{array}{l}{\mathrm{MTBR}}_{(L,N)}={\displaystyle \sum _{n=0}^{\infty }\frac{(n+1)}{(n+2)!}}{\phi }_N{\Delta }_N^{-2}{({\Delta }_{N}L)}^{n+2}{\Delta }_N^0+L{\phi }_N\exp ({\Delta }_{N}L)e\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^N(-\gamma b^{n-1}{G}^{-1}e)\cdot }(1-{\phi }_N\exp ({\Delta }_{N}L)e)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{n=1}^{N-1}(-\gamma b^{n-1}{G}^{-1}e)}\cdot ({\phi }_{n+1}\exp ({\Delta }_{n+1}L)e-{\phi }_n\exp ({\Delta }_{n}L)e)\cdot {\phi }_N\exp ({\Delta }_{N}L)e+(-\delta H^{-1}{e}_{m_1}).\end{array}$$

(20)

4. Optimal Replacement Policy for the System

We know that a renewal reward process is constituted by the successive renewal cycles together with the costs incurred in each cycle. According to the well-known renewal reward theorem (Ross [29]), the long-run average profit per unit time can be represented as

$$C(L,N)=\frac{\mathrm{Expected}\hspace{0.17em}\mathrm{profit}\hspace{0.17em}\mathrm{per}\hspace{0.17em}\mathrm{renewal}\hspace{0.17em}\mathrm{cycle}}{\mathrm{Expected}\hspace{0.17em}\mathrm{length}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{a}\hspace{0.17em}\mathrm{renewal}\hspace{0.17em}\mathrm{cycle}}$$

Now, we include costs, defining the following cost rate:

$c_o$: the operating reward rate while the system is operating;

$c_r$: the repair cost rate while the system is failed;

$c_p$: the replacement cost rate while the system is being replaced;

$B$: the basic cost for replacing the system;

$$\underset{L,N}{\max }C(L,N)=\underset{L,N}{\max }\frac{c_o[E[\min \{{\displaystyle \sum _{n=1}^N{O}_n},L\}]]-c_{r}E[{({\displaystyle \sum _{n=1}^N{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}\le L\}}+{({\displaystyle \sum _{n=1}^{M(L)}{Z}_n})}_{\chi \{{\displaystyle \sum _{n=1}^N{O}_n}>L\}}]-c_{p}E[R]-B}{{\mathrm{MTBR}}_{(L,N)}}.$$

(21)

Substituting Equations (13), (18) and (19) into (21), we obtain

$$\underset{L,N}{\max }C(L,N)=\underset{L,N}{\max }\frac{\left(\begin{array}{l}{c}_o[{\displaystyle \sum _{n=0}^{\infty }\frac{(n+1)}{(n+2)!}}{\phi }_N{\Delta }_N^{-2}{({\Delta }_{N}L)}^{n+2}{\Delta }_N^0+L{\phi }_N\exp ({\Delta }_{N}L)e]\\ -c_r[{\displaystyle \sum _{n=1}^N(-\gamma b^{n-1}{G}^{-1}e)\cdot }(1-{\phi }_N\exp ({\Delta }_{N}L)e)+{\displaystyle \sum _{n=1}^{N-1}(-\gamma b^{n-1}{G}^{-1}e)}\cdot ({\phi }_{n+1}\exp ({\Delta }_{n+1}L)e-{\phi }_n\exp ({\Delta }_{n}L)e)\cdot {\phi }_N\exp ({\Delta }_{N}L)e]\\ -c_p[-\delta H^{-1}{e}_{m_1}]-B\end{array}\right)}{\left(\begin{array}{l}{\displaystyle \sum _{n=0}^{\infty }\frac{(n+1)}{(n+2)!}}{\phi }_N{\Delta }_N^{-2}{({\Delta }_{N}L)}^{n+2}{\Delta }_N^0+L{\phi }_N\exp ({\Delta }_{N}L)e\\ +{\displaystyle \sum _{n=1}^N(-\gamma b^{n-1}{G}^{-1}e)\cdot }(1-{\phi }_N\exp ({\Delta }_{N}L)e)\\ +{\displaystyle \sum _{n=1}^{N-1}(-\gamma b^{n-1}{G}^{-1}e)}\cdot ({\phi }_{n+1}\exp ({\Delta }_{n+1}L)e-{\phi }_n\exp ({\Delta }_{n}L)e)\cdot {\phi }_N\exp ({\Delta }_{N}L)e+(-\delta H^{-1}{e}_{m_1})\end{array}\right)}.$$

(22)

From Equation (22), we can see that it is difficult to develop the optimal solution $(L^{\ast },N^{\ast })$ symbolically due to the high non-linear and complex nature of the optimization problem. Therefore, in the following section, numerical experiments are performed to demonstrate the optimal solution $(L^{\ast },N^{\ast })$ of the long-run average profit rate function $C(L,N)$ that exists and is unique.

5. Numerical Examples

As shown in Figure 2, a micro-engine consists of several orthogonal linear comb drive actuators, which are mechanically joined to a rotating gear. The wear on the rubbing surface between the gear and the pin joint usually causes a broken pin, which is the dominant reason for the failure of micro-engines. Additionally, in shock tests on micro-engines, the fracture of the springs is observed when $k_1$ consecutive external shocks the size of at least $d_1$ or $k_2$ consecutive external shocks the size of at least $d_2$ occur, where $d_1<d_2$, $k_1>k_2$. Therefore, the micro-engine is subject to two competing failure processes: soft failure due to wear degradation and hard failure due to spring fracture caused by external shocks.

The distributions of the random variable involved in the proposed model are given as follows:

(1)

Interarrival time between two successive external shocks: $X_i^n\sim P{H}_c(\beta ,S)$,$i=1,2,\dots ,n=1,2,\dots ,N$,

$$\beta =(1,0),S=\left(\begin{array}{cc}-0.15& 0.15\\ 0& -0.15\end{array}\right),S^0=\left(\begin{array}{c}0\\ 0.15\end{array}\right).$$

(2)

The magnitude of the $i\mathrm{th}$ shock in the $n\mathrm{th}$ sub-cycle: $Y_i^n\sim Exp(5)$$i=1,2,\dots ,n=1,2,\dots ,N$.

(3)

The lifetime of the system due to internal wear degradation in the $n\mathrm{th}$ sub-cycle: $W_n\sim P{H}_c(\alpha ,a^{n-1}T)$, $n=1,2,\dots ,N$,

$$\alpha =(1,0),T=\left(\begin{array}{cc}-0.004& 0.004\\ 0& -0.004\end{array}\right),T^0=\left(\begin{array}{c}0\\ 0.004\end{array}\right),a=1.1.$$

(4)

Repair time after the $n\mathrm{th}$ failure: $Z_n\sim P{H}_c(\gamma ,b^{n-1}G)$, $n=1,2,\dots ,N-1$,

$$\gamma =(1,0),G=\left(\begin{array}{cc}-0.12& 0.12\\ 0& -0.12\end{array}\right),G^0=\left(\begin{array}{c}0\\ 0.12\end{array}\right),b=0.95.$$

(5)

Replacement time: $R\sim P{H}_c(\delta ,H)$,

$$\delta =(1,0),H=\left(\begin{array}{cc}-0.2& 0.2\\ 0& -0.2\end{array}\right),{\mathit{H}}^0=\left(\begin{array}{c}0\\ 0.2\end{array}\right)$$

The other parameters are assumed as: $k_1=3$, $k_2=2$, $d_1=0.5$, $d_2=1$, and the cost parameters are taken as $c_o=210$, $c_r=24$, $c_p=29$, $B=3100$. Figure 3 is the curve representing the profit rate function $C(L,N)$ when $N=10$. As shown in Figure 3, a maximum long-run average profit rate of 934.459 is achieved at ($L=102,N=10$). This means that the system manager should replace the aging system with a new one after the 10th failure or when the total working time of the system exceeds 102 time units, whichever occurs first.

6. Conclusions

In this paper, a multiple competing failure geometric process repair system model under a mixed run shock pattern is proposed and investigated. All random variables involved in the proposed model are phase-type distributed, and the bivariate replacement policy $(L,N)$ is adopted. The density of the PH distributions in the family of the distribution functions defined on the positive real line can be used for approaching any general distributions. We go beyond previous works by not only considering phase-type distributions in all cases and the mixed run shocks, but also adopting the bivariate replacement policy $(L,N)$. The long-run average profit rate for the system is obtained by employing the closure property of the PH distribution. It is worth noting that the expression of the long-run average profit rate for the system $C(L,N)$ is high non-linear and complex. Therefore, we never obtained the optimal analytical solution of the corresponding optimization problem. The dimension of the constructed system model will increase sharply due to the use of PH distributions. To find a way to reduce the dimension is also an important issue for future study. In addition, other mixed shock models under MDCFP can be studied. For example, for two fixed critical values $\delta$ and $L$, the system fails when k out of interarrival times between two successive shocks are less than $\delta$, or the cumulative magnitude of the shock is larger than $L$.