2. Preliminaries and model description

We begin by introducing the notation that will be used throughout the paper. For any underlying random variable U, we denote the probability density function by $f_{U}(\cdot)$ , the cumulative distribution function by $F_{U}(\cdot)$ , the reliability function by $\bar{F}_{U}(\cdot)$ , the failure rate function by $r_{U}(\cdot)$ , the Laplace transform function by $L_{U}(\cdot)$ , and the mean by $\mathbb{E}(U)$ . Further, we denote the set of natural numbers by $\mathbb{N}$ and the set of real numbers by $\mathbb{R}$ . The random variable $X_{i}$ represents the time between the $(i-1)$ th and ith shocks, and $Y_{i}$ denotes the magnitude (or damage) caused by the ith shock, $i \in \mathbb{N}$ . Further, the random variables T and N represent the lifetime of the system and the number of shocks until system failure, respectively. We denote the derivative of a function g(x) with respect to x by g ^′(x), the absolute value of a real number z by $|z|$ , and the identity matrix by $\boldsymbol{I}$ . The Kronecker product of any two matrices $\pmb{E}$ and $\pmb{F}$ is denoted by $\pmb{E} \otimes \pmb{F} = [E_{ij}\pmb{F}]$ , and $\mathrm{diag}(x_{1},x_{2}, \ldots,x_{d})$ denotes the diagonal matrix with ith diagonal entry equal to $x_{i}$ , $i = 1,2,\ldots,d$ .

2.1. Preliminaries

For the stochastic objects considered here, we present the relevant definitions in this subsection.

Definition 1. For two matrices $\pmb{S}_{1}$ and $\pmb{S}_{2}$ , we write $\pmb{S}_{1} \geq \pmb{S}_{2}$ (respectively $\pmb{S}_{1} > \pmb{S}_{2}$ ) if $\pmb{S}_{1}$ and $\pmb{S}_{2}$ are real and every entry of $\pmb{S}_{1} - \pmb{S}_{2}$ is non-negative (respectively positive).

Definition 2. A real matrix $\pmb{S}$ is said to be monotone if $\pmb{S} \pmb{x} \geq 0$ implies $\pmb{x} \geq 0$ .

Definition 3. A matrix $\pmb{S}$ is called an M-matrix if all off-diagonal elements are less than or equal to zero, and all eigenvalues have positive real parts [Reference He17, p. 22].

Remark 1. Any non-singular M-matrix is a monotone matrix, and its inverse is a non-negative matrix [Reference Johnson, Smith and Tsatsomeros18].

Definition 4. A non-negative random variable X is said to have a phase-type distribution if

\begin{equation*} F_X(x) = 1 - {\alpha}{\mathrm{e}}^{\pmb{S}x}\boldsymbol{e} = 1 - {\alpha}\Bigg(\sum_{i=0}^{\infty}\frac{x^{i}}{i!}\pmb{S}^{i}\Bigg)\boldsymbol{e}, \quad x\geq 0, \end{equation*}

where:

• • $\boldsymbol{e}$ is the column vector of order m, where $m \in \mathbb{N}$ , with all entries equal to one;

• • ${\alpha}$ is a row vector of order m with all elements non-negative, and ${\alpha}\boldsymbol{e} = 1$ ; and

• • $\pmb{S}$ is a matrix of order m such that all the entries are non-negative except the diagonal elements, which are negative; all the row sums are less than or equal to zero; and $\pmb{S}$ is non-singular [Reference He17, p. 10].

We call $\pmb{S}$ and the pair $({\alpha},\pmb{S})$ the phase-type generator and the phase-type representation of order m, respectively. The notation $X \sim \mathrm{PH}_{m}({\alpha},\pmb{S})$ means that X follows the phase-type distribution with the representation $({\alpha},\pmb{S})$ of order m. The expectation of X can be derived from

(1) \begin{equation} \mathbb{E}(X) = -({\alpha} \pmb{S}^{-1} \boldsymbol{e}).\end{equation}

The class of phase-type distributions includes many well-known probability distributions. For example, a phase-type distribution with the representation $({\alpha} = 1, \pmb{S} = [\!-\!\lambda])$ corresponds to the exponential distribution with parameter $\lambda$ . Additionally, the Erlang distribution, generalized Erlang distribution, and Coxian distribution are notable special cases of phase-type distributions; for the phase-type representation of these distributions, see [Reference He17, Example 1.2.3].

Remark 2. Let $\pmb{S}$ be a phase-type generator. Then all the eigenvalues of $\pmb{S}$ have negative real part [Reference Bladt and Nielsen2]. Therefore, from Definitions 3 and 4, for $t \geq 0$ , $(t\boldsymbol{I} - \pmb{S})$ is an M-matrix. Hence, from Remark 1, it is monotone.

Definition 5. A discrete phase-type distribution can be viewed as the distribution of the time to absorption in an absorbing Markov chain. Let $\boldsymbol{M}$ be a square matrix of size m, $m \in \mathbb{N}$ , that includes the transition probabilities among the m transient states, and $\pmb{\mu} = (\mu_{1},\mu_{2},\ldots,\mu_{m})$ be a non-negative vector of size m such that $\sum_{i=1}^{m}\mu_{i} = 1$ . Further, let the matrix $\boldsymbol{I} - \boldsymbol{M}$ be non-singular. Then a random variable N is said to have the discrete phase-type distribution, denoted by $N \sim \mathrm{DPH}_{m}(\pmb{\mu}, \boldsymbol{M})$ , if

(2) \begin{equation} \mathbb{P}(N = n) = \pmb{\mu}\boldsymbol{M}^{n-1}(\boldsymbol{I} - \boldsymbol{M})\pmb{e}, \quad n \in \mathbb{N}, \end{equation}

where $\pmb{e}$ is the column vector of size m with all elements being one [Reference Goyal, Hazra and Finkelstein14].

Definition 6. Given matrices $\{\boldsymbol{D}_{0},\boldsymbol{D}_{1}\}$ of order m, where $m \in \mathbb{N}$ , all the elements of the two matrices are non-negative except the diagonal elements of $\boldsymbol{D}_{0}$ , which are negative, and $\boldsymbol{D} = \boldsymbol{D}_{0} + \boldsymbol{D}_{1}$ is an irreducible infinitesimal generator of a Markov process, consider an m-state Markov renewal process $\{(I_{k}, X_{k}), k = 0,1,2,\ldots\}$ . The random variables $I_{k}$ and $X_{k}$ are the states and the sojourn times in the states, respectively. The transition probability matrix with entries (i, j) is given by

$$ \mathbb{P}(I_{k} = j, X_{k} \leq t\mid I_{k-1} = i) = \bigg(\int_{0}^{t}{\mathrm{e}}^{\boldsymbol{D}_{0}x}\boldsymbol{D}_{1}\,{\mathrm{d}} x\bigg)_{ij}, \quad t \geq 0. $$

The Markovian arrival process is defined as the Markov renewal process and a transition probability matrix of the stated particular form. We call the pair $(\boldsymbol{D}_{0},\boldsymbol{D}_{1})$ the representation of the Markovian arrival process. If $\pmb{\pi}$ denotes the initial vector, then we denote the Markovian arrival process by $\mathrm{MAP}_{m}(\pmb{\pi},\boldsymbol{D}_{0}, \boldsymbol{D}_{1})$ [Reference He17, p. 111].

The Markovian arrival process generalizes the Poisson process, the Markov modulated Poisson process, and the phase-type renewal process [Reference Latouche, Remiche and Taylor21]. The joint distribution of $(X_{1},X_{2},\ldots,X_{k})$ for a $\mathrm{MAP}_{m}(\pmb{\pi}, \boldsymbol{D}_{0}, \boldsymbol{D}_{1})$ process is given by

(3) \begin{equation} f_{(X_{1},X_{2},\ldots,X_{k})}(x_{1},x_{2},\ldots,x_{k}) = \pmb{\pi}{\mathrm{e}}^{\boldsymbol{D}_{0}x_{1}}\boldsymbol{D}_{1} {\mathrm{e}}^{\boldsymbol{D}_{0}x_{2}}\cdots\boldsymbol{D}_{1}{\mathrm{e}}^{\boldsymbol{D}_{0}x_{k}}\boldsymbol{D}_{1}\boldsymbol{e}\end{equation}

for $x_{i} \geq 0$ , $i = 1,2,\ldots,k$ . Clearly, all the inter-arrival times $X_{i}$ are dependent, non-identical, and follow a phase-type distribution with generator $\boldsymbol{D}_{0}$ [Reference Bladt and Nielsen2, p. 526].

2.2. Model description

Consider a system subject to shocks over time which are the only factor contributing to its failure. Let $X_{i}$ denote the time between the $(i-1)$ th and ith shocks, and $Y_{i}$ represent the magnitude of the ith shock, $i \in \mathbb{N}$ . Let the random variable N indicate the number of shocks until the system fails. Then the system’s lifetime can be represented by the random sum

(4) \begin{equation} T = \sum_{i=1}^{N} X_{i}.\end{equation}

Note that this is valid for any system working in a random environment. However, by assuming an arbitrary distribution of N, we cannot proceed much further and derive in full generality the distribution of the lifetime of a system vulnerable to shocks. We assume that $N \sim \mathrm{DPH}_{m}(\pmb{\mu},\boldsymbol{M})$ . The discrete phase-type distribution is reasonably general since some popular waiting time distributions, e.g. geometric, geometric distribution of order k, and negative binomial, are special cases. Consequently, many shock models have a discrete phase-type distribution of N (see Section 4). Furthermore, this distribution is dense in the class of distributions in $\mathbb{N}$ [Reference Latouche and Ramaswami20, p. 54]. Therefore, it is a somewhat broad assumption.

As with the distribution of N, we cannot derive the system lifetime by considering an arbitrary shock process. We assume that shocks occur on the system according to the Markovian arrival process, which is dense in the class of point processes on the real line [Reference Bladt and Nielsen2, p. 708] and provides great flexibility in the model. Further, we assume that $X_{i}$ and N are independent, $i \in \mathbb{N}$ . Lastly, we assume that the $Y_{i}$ are independent and identically distributed random variables, and independent of the $X_{i}$ .

We can observe that the class we have defined is very broad but not complete (i.e. it does not contains all shock models, e.g. the $\delta$ -shock model). The $\delta$ -shock model states that the system fails when the time lag between two successive shocks is less than a predetermined threshold $\delta$ [Reference Goyal, Hazra and Finkelstein13]. For this model, $\{N = n\}$ if and only if $\{X_{1} > \delta, X_{2} > \delta, \ldots, X_{n-1} > \delta, X_{n} \leq \delta\}$ . Since N and the $X_{i}$ are dependent, the $\delta$ -shock model is not part of the proposed class.

3. Reliability and mean residual lifetime of the system

We begin this section with the following lemma which is used to prove the main result of this section.

In what follows, we obtain the distribution of the system’s lifetime for the proposed class of shock models.

Proof. Clearly, from (4), and by using the assumption of independence between the $X_{i}$ and N, we can write

(5) \begin{equation} L_{T}(t) = \mathbb{E}({\mathrm{e}}^{-T t}) = \sum_{n=1}^{\infty}\mathbb{E}\big({\mathrm{e}}^{-\left(\sum_{i=1}^{n}X_{i}\right)t}\big) \mathbb{P}(N = n). \end{equation}

Now, $\mathbb{E}\big({\mathrm{e}}^{-\left(\sum_{i=1}^{n}X_{i}\right)t}\big) = \int_{0}^{\infty}\cdots\int_{0}^{\infty}{\mathrm{e}}^{\left(-\sum_{i=1}^{n}x_{i}\right)t}f_{X_{1},\ldots,X_{n}}(x_{1},\ldots,x_{n})\, {\mathrm{d}} x_{1} \cdots {\mathrm{d}} x_{n}$ . From (3), we can write

\begin{align*} \mathbb{E}\big({\mathrm{e}}^{-\left(\sum_{i=1}^{n}X_{i}\right)t}\big) & = \int_{0}^{\infty}\cdots\int_{0}^{\infty}{\mathrm{e}}^{\left(-\sum_{i=1}^{n}x_{i}\right)t} \big(\pmb{\pi}{\mathrm{e}}^{\boldsymbol{D}_{0}x_{1}}\boldsymbol{D}_{1}{\mathrm{e}}^{\boldsymbol{D}_{0}x_{2}}\cdots\boldsymbol{D}_{1} {\mathrm{e}}^{\boldsymbol{D}_{0}x_{n}}\boldsymbol{D}_{1}\boldsymbol{e}\big)\,{\mathrm{d}} x_{1}\cdots{\mathrm{d}} x_{n} \\ & = \int_{0}^{\infty}\cdots\int_{0}^{\infty}\pmb{\pi}{\mathrm{e}}^{-x_{1}\left(t\boldsymbol{I}-\boldsymbol{D}_{0}\right)} \boldsymbol{D}_{1}{\mathrm{e}}^{-x_{2}\left(t\boldsymbol{I}-\boldsymbol{D}_{0}\right)}\cdots\boldsymbol{D}_{1} {\mathrm{e}}^{-x_{1}\left(t\boldsymbol{I}-\boldsymbol{D}_{0}\right)}\boldsymbol{D}_{1}\boldsymbol{e}\,{\mathrm{d}} x_{1}\cdots{\mathrm{d}} x_{n} \\ & = \pmb{\pi}(t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1}(t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1} \cdots\boldsymbol{D}_{1}(t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1}\boldsymbol{e} \\ & = \pmb{\pi}((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1})^{n}\boldsymbol{e}. \end{align*}

Now, by using this equality and (2) in (5), we get

\begin{align*} L_{T}(t) & = \sum_{n=1}^{\infty}(\pmb{\pi}((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1})^{n}\boldsymbol{e}) (\pmb{\mu}\boldsymbol{M}^{n-1}(\boldsymbol{I}-\boldsymbol{M})\boldsymbol{e}) \\ & = \sum_{n=1}^{\infty}(\pmb{\pi}\otimes\pmb{\mu})\{((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1} \boldsymbol{D}_{1})^{n-1}\otimes\boldsymbol{M}^{n-1}\}((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1} \boldsymbol{D}_{1}\boldsymbol{e}\otimes(\boldsymbol{I}-\boldsymbol{M})\boldsymbol{e}) \\ & = (\pmb{\pi}\otimes\pmb{\mu})\Bigg[\sum_{n=1}^{\infty}((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1} \boldsymbol{D}_{1}\otimes\boldsymbol{M})^{n-1}\Bigg]((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\otimes \boldsymbol{I})(\boldsymbol{D}_{1}\boldsymbol{e}\otimes(\boldsymbol{I}-\boldsymbol{M})\boldsymbol{e}). \end{align*}

The series term in this equation converges to $(\boldsymbol{I}-((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1} \otimes\boldsymbol{M}))^{-1}$ if and only if all the eigenvalues of the matrix $(t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1}\otimes\boldsymbol{M}$ has modulus strictly less than 1. As all the eigenvalues for the matrix $\boldsymbol{M}$ of a discrete phase-type distribution have modulus strictly less than one [Reference Bladt and Nielsen2, p. 32], from Lemma 1, all the eigenvalues of the matrix $(t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1}\boldsymbol{D}_{1}\otimes\boldsymbol{M}$ lie strictly within the unit circle. Hence, by using $\boldsymbol{D}_{1}\boldsymbol{e} = -\boldsymbol{D}_{0}\boldsymbol{e}$ , we get

\begin{align*} L_{T}(t) & = (\pmb{\pi}\otimes\pmb{\mu})(\boldsymbol{I}-((t\boldsymbol{I}-\boldsymbol{D}_{0})^{-1} \boldsymbol{D}_{1}\otimes\boldsymbol{M}))^{-1}((t\boldsymbol{I}-\boldsymbol{D}_{0})\otimes\boldsymbol{I})^{-1} (\boldsymbol{D}_{1}\boldsymbol{e}\otimes(\boldsymbol{I}-\boldsymbol{M})\boldsymbol{e}) \\ & = (\pmb{\pi}\otimes\pmb{\mu})(t\boldsymbol{I}-(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M}))^{-1}(\boldsymbol{D}_{1}\boldsymbol{e}\otimes\boldsymbol{e} - \boldsymbol{D}_{1}\boldsymbol{e}\otimes\boldsymbol{M}\boldsymbol{e}) \\ & = (\pmb{\pi}\otimes\pmb{\mu})(t\boldsymbol{I}-(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M}))^{-1}(\!-\!\boldsymbol{D}_{0}\boldsymbol{e}\otimes\boldsymbol{e} - \boldsymbol{D}_{1}\boldsymbol{e}\otimes\boldsymbol{M}\boldsymbol{e}) \\ & = (\pmb{\pi}\otimes\pmb{\mu})(t\boldsymbol{I}-(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M}))^{-1}(\!-\!\boldsymbol{D}_{0}\otimes\boldsymbol{I} - \boldsymbol{D}_{1}\otimes\boldsymbol{M})\boldsymbol{e} \\ & = (\pmb{\pi}\otimes\pmb{\mu})(t\boldsymbol{I}-(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M}))^{-1}(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M})^{0}, \end{align*}

where $(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M})^{0} = -(\boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M})\boldsymbol{e}$ . Which is the Laplace transform of a phase-type distribution with the representation $(\pmb{\pi}\otimes\pmb{\mu},\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})$ [Reference He17, p. 18]. Consequently, $T \sim \mathrm{PH}(\pmb{\pi}\otimes\pmb{\mu}, \boldsymbol{D}_{0}\otimes\boldsymbol{I} + \boldsymbol{D}_{1}\otimes\boldsymbol{M})$ .

Theorem 1 shows that the system’s lifetime follows the well-known phase-type distribution for the defined class of shock models. Many shock models and their applications have been studied widely in the literature for phase-type renewal processes (see [Reference Eryilmaz and Kan9–Reference Gong, Xie and Yang11, Reference Lyu, Qu, Ma, Wang and Yang23, Reference Ozkut30, Reference Wang, Zhao, Wu and Wang36, Reference Zhao, Wang, Wang and Fan38], to name a few). We believe the theorem’s conclusion can be applied to all of these current and upcoming studies, making this conclusion of considerable importance.

Since, for the defined class of shock models, the system’s lifetime follows the phase-type distribution, from Definition 4, the system’s reliability function is given by

(6) \begin{equation} \bar{F}_{T}(t) = (\pmb{\pi}\otimes\pmb{\mu}){\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})t} \boldsymbol{e},\end{equation}

and the failure rate is given by

(7) \begin{equation} r_{T}(t) = \frac{-(\pmb{\pi}\otimes\pmb{\mu}){\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})t} (\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})\boldsymbol{e}} {(\pmb{\pi}\otimes\pmb{\mu}){\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})t} \boldsymbol{e}}.\end{equation}

For any phase-type distribution $Y \sim \mathrm{PH}_{m}({\alpha}, \pmb{S})$ ,

$$(Y - y\mid Y > y) \sim \mathrm{PH}\bigg(\frac{{\alpha}{\mathrm{e}}^{\pmb{S}y}}{{\alpha}{\mathrm{e}}^{\pmb{S}y}\boldsymbol{e}},\pmb{S}\bigg)$$

(see [Reference He17]). Thus, the system’s mean remaining lifetime is given by

(8) \begin{equation} \mathbb{E}(T - t\mid T > t) = -\bigg(\frac{(\pmb{\pi}\otimes\pmb{\mu}){\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})t}} {(\pmb{\pi}\otimes\pmb{\mu}){\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})t} \boldsymbol{e}}\bigg) (\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})^{-1}\boldsymbol{e}.\end{equation}

4. Optimal age replacement policy via phase-type distribution

The classical age replacement policy states that the system is replaced when it fails or reaches age x, whichever comes first. Let $c_{1}$ and $c_{2}$ represent, respectively, the costs of replacing working and failed systems. We assume that $c_{1} < c_{2}$ since a failure results in an additional penalty. Then, as a function of the replacement age x, the average cost per unit of time is given by

(9) \begin{equation} C(x) = \frac{c_{1}\mathbb{P}(T > x) + c_{2}\mathbb{P}(T \leq x)}{\mathbb{E}(\mathrm{min}(T,x))} = \frac{c_{2} + (c_{1} - c_{2}) P(T > x)}{\mathbb{E}(\mathrm{min}(T,x))}.\end{equation}

To evaluate the function C(x) we need expressions for $\mathbb{P}(T > x)$ and $\mathbb{E}(min(T,x))$ . Note that, for a phase-type distribution Y which has a phase-type representation $({\alpha}, \pmb{S})$ , i.e. $Y \sim \mathrm{PH}_{m}({\alpha}, \pmb{S})$ , $\mathbb{E}(\min(Y,x)) = {\alpha}\pmb{S}^{-1}{\mathrm{e}}^{\pmb{S}x}\boldsymbol{e} - {\alpha}\pmb{S}^{-1}\boldsymbol{e}$ . Hence, from this equation and Theorem 1,

\begin{align*}\mathbb{E}(\!\min(T,x)) &= (\pmb{\pi}\otimes\pmb{\mu})(\boldsymbol{D}_{0}\otimes\boldsymbol{I} +\boldsymbol{D}_{1}\otimes\boldsymbol{M})^{-1}{\mathrm{e}}^{(\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})x}\pmb{e}\\& - (\pmb{\pi}\otimes\pmb{\mu})(\boldsymbol{D}_{0}\otimes\boldsymbol{I} +\boldsymbol{D}_{1}\otimes\boldsymbol{M})^{-1}\pmb{e}\end{align*}

for $x \geq 0$ . Now, by using this and (6) we can obtain value of the cost function C(x) numerically.

We get a phase-type distributed system lifetime for the defined class of shock models. This leads us to consider the lifetime random variable T to be phase-type distributed in the age replacement policy. This policy was studied in [Reference Eryilmaz8] by considering the system lifetime distribution to be a discrete phase-type distribution, while we are considering a continuous phase-type distribution.

For a given shock model and costs $c_{1}$ and $c_{2}$ , the problem is finding the point $x^{*}$ that minimizes the function C(x). For the classical model (9), the optimal solution exists for the increasing failure rate with the condition that it tends to a constant. Our model also agrees with this, which we prove in the following theorem.

Theorem 2. Assume that the lifetime of the system follows a phase-type distribution, i.e. $T \sim \mathrm{PH}_{m}({\alpha},\pmb{S})$ . Let $\theta$ be the dominating eigenvalue of the matrix $\pmb{S}$ . Let $\sigma = \mathbb{E}(T)$ . Assume that $r_{T}(\cdot)$ is strictly increasing and continuous. If $|\theta| > c_{2}/\sigma(c_{2} - c_{1})$ then there exists a finite and unique $x^{*}$ ( $0 < x^{*} < \infty$ ) that satisfies

$$r_{T}(x) \int_{0}^{x} \bar{F}_{T}(t)\,{\mathrm{d}} t - F(x) = \frac{c_{1}}{c_{2} - c_{1}},$$

and the resulting expected cost rate is $C(x^{*}) = (c_{2} - c_{1}) r_{T}(x^{*})$ . Further, if $|\theta| \leq c_{2}/\sigma(c_{2} - c_{1})$ then $x^{*} = \infty$ , i.e. a unit is replaced only at failure, and the expected cost rate is given in (9).

Proof. First, note that, for a phase-type distribution, the dominating eigenvalue is always real and negative [Reference Bladt and Nielsen2, p. 203]. From [Reference Nakagawa28, Theorem 3.1], a sufficient condition that $\lim_{x \rightarrow \infty} C(x) > C(x)$ for some x is that $\lim_{t \rightarrow \infty} r_{T}(t) > c_{2}/\sigma(c_{2} - c_{1})$ . Note that the survival function of a phase-type distribution T can be written as

(10) \begin{equation} \bar{F}_{T}(t) = \sum_{k=1}^{l}\sum_{r=0}^{m_{k} - 1}s_{k_{r}}\frac{t^{r}}{r!}{\mathrm{e}}^{\theta_{k}t}, \end{equation}

where the $\theta_{i}$ are the eigenvalues of $\pmb{S}$ with multiplicities $m_{1}, m_{2}, \ldots, m_{l}$ such that $\sum_{i=1}^{l} m_{i} = m$ , while the $s_{k_{r}}$ are constants [Reference Bladt and Nielsen2]. From (10), we can see that the failure rate function of T is always asymptotically constant and is $|\theta|$ , i.e. $\lim_{t \rightarrow \infty}r_{T}(t) = |\theta|$ . Then our result follows directly from [Reference Nakagawa28, Theorem 3.2].

Since the failure rate function, $r_{T}(\cdot)$ in (7), is matrix based, its monotonic behavior can be seen by plotting it in any mathematical software. However, in the following two theorems, adequate criteria for the monotonicity of the failure rate of a phase-type distribution are derived for two exceptional instances.

Proof. Define $\beta_{i} = {\alpha}\boldsymbol{C}_{i}\pmb{R}_{i}\boldsymbol{e}$ , $i = 1,2$ . Since the dominating eigenvalue of a phase-type distribution is always real, this implies that both the eigenvalues will be real. Which means $\pmb{P}$ and $\pmb{P}^{-1}$ will be real matrices, as, consequently, are the $\pmb{C}_{i}$ and $\pmb{R}_{i}$ . Therefore, $\beta_{i} \in \mathbb{R}$ , $i = 1,2$ . Now, consider

\begin{align*} \bar{F}_{T}(t) = {\alpha}{\mathrm{e}}^{\pmb{S}t}\boldsymbol{e} = {\alpha}\pmb{P}\,{\mathrm{diag}}({\mathrm{e}}^{\theta_{1}t},{\mathrm{e}}^{\theta_{2}t})\pmb{P^{-1}}\boldsymbol{e} & = {\alpha}(\boldsymbol{C}_{1},\boldsymbol{C}_{2})\,\mathrm{diag}({\mathrm{e}}^{\theta_{1}t},{\mathrm{e}}^{\theta_{2}t}) \begin{pmatrix} \pmb{R}_{1} \\ \pmb{R}_{2} \end{pmatrix} \boldsymbol{e} \\ & = ({\alpha}\pmb{C}_{1},{\alpha}\boldsymbol{C}_{2})\, \mathrm{diag}({\mathrm{e}}^{\theta_{1}t},{\mathrm{e}}^{\theta_{2}t}) \begin{pmatrix} \pmb{R}_{1}\boldsymbol{e} \\ \pmb{R}_{2}\boldsymbol{e} \end{pmatrix} \\ & = {\mathrm{e}}^{\theta_{1}t}({\alpha}\boldsymbol{C}_{1}\pmb{R}_{1}\boldsymbol{e}) + {\mathrm{e}}^{\theta_{2}t}({\alpha}\boldsymbol{C}_{2}\pmb{R}_{2}\boldsymbol{e}) = \beta_{1}{\mathrm{e}}^{\theta_{1}t} + \beta_{2}{\mathrm{e}}^{\theta_{2}t}. \end{align*}

Thus, $f_{T}(t) = -\beta_{1}\theta_{1}{\mathrm{e}}^{\theta_{1}t} - \beta_{2}\theta_{2}{\mathrm{e}}^{\theta_{2}t}$ . Simple calculation gives $r^{\prime}_{T}(t) = -\beta_{1}\beta_{2}(\theta_{1} - \theta_{2})^{2}{\mathrm{e}}^{(\theta_{1} + \theta_{2})t}$ , which is positive if and only if $\beta_{1} \beta_{2} < 0$ . Now, observe that

$$ \beta_{1} + \beta_{2} = {\alpha}\boldsymbol{C}_{1}\pmb{R}_{1}\boldsymbol{e} + {\alpha}\boldsymbol{C}_{2}\pmb{R}_{2}\boldsymbol{e} = {\alpha}(\boldsymbol{C}_{1},\boldsymbol{C}_{2}) \begin{pmatrix} \pmb{R}_{1}\\ \pmb{R}_{2} \end{pmatrix}\boldsymbol{e} = {\alpha}\pmb{P}\pmb{P^{-1}}\boldsymbol{e} = 1. $$

The result follows using this equality in $\beta_{1}\beta_{2} < 0$ .

5. Special shock models

In this section, we briefly review models that fall under the defined category of shock models for which N follows a discrete phase-type distribution. Many such models have been studied in the literature, more than we can list here. Therefore, we illustrate some representative shock models. Further, we also show that the extreme shock model with a finite number of shocks belongs to the defined class of shock models (which has not been observed in the literature).

Since all the shock models discussed in this section fall within the defined class, the system’s lifetime distribution and other reliability measures can be derived using Theorem 1 and (6)–(8).

5.1. Extreme shock model

One of the shock models that has drawn considerable attention in the reliability literature is the extreme shock model. Here, the system fails when an individual shock exceeds some given level, e.g. when a crack is too deep, it might cause a vehicle axle to fail. That is, if the magnitude of the ith shock exceeds a predetermined threshold $\gamma$ ( $Y_{i} > \gamma$ ) the system fails at the ith shock, $i \in \mathbb{N}$ [Reference Gut and Hüsler16]. In this case,

\begin{equation*} \mathbb{P}(N = n) = \mathbb{P}(Y_{1} \leq \gamma, Y_{2} \leq \gamma, \ldots, Y_{n-1} \leq \gamma, Y_{n} > \gamma)= p(1-p)^{n-1},\end{equation*}

where $p = \mathbb{P}(Y_{1} > \gamma)$ . Clearly, $N \sim \mathrm{DPH}_{1}(1,1 - p)$ . Other than the classical extreme shock model, some generalized extreme shock models with discrete phase-type representation of N have also been put forward in the literature [Reference Bozbulut and Eryilmaz3].

5.2. A cumulative shock model

A system fails in the cumulative shock model if the total damage from shocks reaches a predetermined threshold value $\gamma$ . This model has been studied in a lot of real applications. For instance, it is possible to imagine each electrical device connected to an electric system as a specific type of shock; the power of the device is comparable to the magnitude of shocks. The system will be overloaded, resulting in a power failure, when the total shock size surpasses a certain level. Let the random variable N(t) denote the number of shocks by time t, $t \geq 0$ . Suppose that the ith shock increases the damage of the system by independent and identically distributed random increments $Y_{i}$ , $i \in \mathbb{N}$ , which are independent of the shock process $\{N(t)\colon t \geq 0\}$ . Therefore, the system’s total damage at time t is determined by the random sum $Y(t) = \sum_{i=1}^{N(t)} Y_{i}$ , where $Y(t) = 0$ when $N(t) = 0$ . Suppose a system fails at time t when the cumulative damage surpasses a random boundary B, $Y(t) > B$ (note that in this model, the threshold value $\gamma$ is considered random). Therefore, in this case, $\bar{F}_{T}(t) = \mathbb{P}(Y(t) \leq B)$ . Let B follow an exponential distribution, and be independent of the shock process and $Y_{i}$ , $i \in \mathbb{N}$ . Then, due to the memoryless property of the exponential distribution, a shock that occurs at time t causes the system to fail with probability $\mathbb{P}(B < Y_{N(t)}) = \mathbb{P}(B < Y_{1})$ [Reference Cha and Finkelstein4]. Thus, in this setting,

\begin{equation*} \mathbb{P}(N = n) = \mathbb{P}(B < Y_{1})\mathbb{P}(B \geq Y_{1})^{n-1} \sim \mathrm{DPH}_{1}(1, \mathbb{P}(B \geq Y_{1})).\end{equation*}

A cumulative shock model studied in [Reference Lyu, Qu, Ma, Wang and Yang23] is also a special case of the defined class of shock models.

5.3. Run shock model

According to the run shock model, a system collapses when the magnitudes of k successive shocks go over a specific threshold $\gamma$ [Reference Mallor and Omey24]. In this case, N can be represented as a discrete phase-type distribution i.e. $N \sim \mathrm{DPH}_{k}(\pmb{\mu},\boldsymbol{M})$ , where $\pmb{\mu} = (1,0,\ldots,0)$ and

$$ \boldsymbol{M} =\left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1 - p & p&0 & \cdots & 0 \\1 - p & 0 & p & \cdots & 0 \\\vdots& \vdots & \vdots & \ddots & \vdots \\1-p & 0 & 0 & \cdots & p\\1 - p&0 & 0 & \cdots & 0\end{array}\right)$$

(see [Reference Tank and Eryilmaz35]). Many generalizations of the run shock model have appeared in the literature. For example, [Reference Eryilmaz and Tekin10] studied a shock model by combining the classical run and extreme shock models, according to which a system fails if either k consecutive shocks occur of size at least $\gamma_{1}$ or a single shock occurs of size at least $\gamma_{2}$ ( $\gamma_{1} < \gamma_{2}$ ). Such a mixed shock model helps model the failure pattern of various devices or systems. For example, a single mechanical stress of small size may not cause the failure of a device, while a certain number of consecutive stresses of small sizes or a single large stress may cause the failure of the device [Reference Eryilmaz and Tekin10]. Another generalization was studied in [Reference Gong, Xie and Yang11], which considered two critical levels $\gamma_1$ and $\gamma_2$ such that $\gamma_{1} < \gamma_{2}$ . According to this model, the system fails if either at least $k_{1}$ consecutive shocks with magnitude above $\gamma_{1}$ or $k_{2}$ consecutive shocks with magnitude above $\gamma_{2}$ occur. This model can be helpful in the healthcare management area, where prostate-specific antigen levels in the treatment of patients was divided into three levels with two thresholds; the type of situation determined the choice of therapy [Reference Gong, Xie and Yang11]. A generalized run shock model was studied in [Reference Yalcin, Eryilmaz and Bozbulut37] by considering system failure due to r non-overlapping runs of k consecutive critical shocks, which helps describe the performance of r cold-standby identical systems. In all the aforementioned shock models, N follows the discrete phase-type distribution.

5.4. Extreme shock model with a change point

This model, in which the system fails according to the extreme shock model with threshold value $\gamma$ , was studied in [Reference Eryilmaz and Kan9]. However, after a random number of shocks, the magnitude distribution changes. A model like this is helpful in real-world applications as a sudden change in environmental factors could result in a greater shock. Consider wind turbine blades that are subject to external shocks; a huge shock may cause a blade crack, which stops the turbine. In fact, the wind pattern or temperature fluctuates over a year, and in these environmental conditions a sudden change could happen. Statistically, the change point is random, and the probabilistic laws of shock sizes might be significantly different before and after the change point. Thus, this shock model is able to simulate such real-world scenarios. Let K represent the number of shocks until the change point. The distribution of the magnitude changes after a random number of shocks. That is, $Y_{1}, \ldots, Y_{K}$ have a common distribution $F_{1}$ such that the shock at K or before is below the threshold $\gamma$ with probability $p_{1} = \mathbb{P}(Y_{i} < \gamma) = F_{1}(\gamma)$ , $i = 1,2,\ldots, K$ , and the shock sizes $Y_{K+1}, Y_{K+2}, \ldots$ after the Kth shock have another common distribution $F_{2}$ such that $p_{2} = \mathbb{P}(Y_{i} < \gamma) = F_{2}(\gamma)$ , $i = K+1, K+2, \ldots$

This configuration has two stages in terms of the shock size that the system experiences. In the first stage, the sizes of the shocks may be relatively small; in the second stage, the system can encounter more severe shocks due to a sudden change in environmental conditions [Reference Eryilmaz and Kan9]. In most cases, this change point is random due to uncertainty. In this case, we have $p_{1} > p_{2}$ . Assume that K follows a geometric distribution with probability mass function $\mathbb{P}(K = k) = q(1-q)^{k-1}$ , $k \in \mathbb{N}$ . In this case, [Reference Eryilmaz and Kan9] showed that, if $p_{1}(1-q) \neq p_{2}$ , then

$$\mathbb{P}(N = n) = \varepsilon[p_{1}(1-q)]^{n-1}[1 - p_{1}(1-q)] + (1 - \varepsilon)p_{2}^{n-1}(1 - p_{2}),$$

where

$$\varepsilon = \frac{p_{1} - p_{2}}{p_{1}(1 - q) - p_{2}}, \qquad p_{1}(1 - q) \neq 1, \qquad p_{1}(1- q) \neq p_{2}.$$

Note that N is a mixture of two geometric distributions with mean $1/(1 - p_{1}(1 - q))$ and $1/(1 - p_{2})$ , which means N follows a discrete phase-type distribution. In this case, $N \sim \mathrm{DPH}_{2}(\pmb{\mu},\boldsymbol{M})$ , where $\pmb{\mu} = (\varepsilon,1-\varepsilon)$ and

$$\boldsymbol{M} = \left(\begin{array} {c@{\quad}c} p_{1}(1-q) & 0 \\ 0 & p_{2}\end{array}\right).$$

Further, if $p_{1}(1-q) = p_{2}$ , then $\mathbb{P}(N = n) = \varepsilon (n-1) p_{2}^{n-2} (1 - p_{2})^{2} + (1 - \varepsilon) p_{2}^{n-1} (1 - p_{2})$ , $n \geq 1$ , where $\varepsilon = ({p_{1} - p_{2}})/({1 - p_{2}})$ . In this case, N is a mixture of geometric distribution and negative binomial, which means N follows a discrete phase-type distribution. In this case, $N \sim \mathrm{DPH}_{3}(\pmb{\mu},\boldsymbol{M})$ , where $\pmb{\mu} = (\varepsilon,0,1-\varepsilon)$ and

$$\boldsymbol{M} = \left(\begin{array} {c@{\quad}c@{\quad}c} p_{2} & 1 - p_{2} & 0 \\ 0 & p_{2} & 0 \\ 0 & 0 & p_{2}\end{array}\right).$$

5.5. Extreme shock model with finite number of shocks

We assume that the system can bear at most J shocks. In this model, the system fails according to the extreme shock model; if the system is not failed by the first J shocks, then it will fail on the arrival of the $(J+1)$ th shock. Let $p_{i}$ be the probability that the system fails at the ith shock, $i = 1,2, \ldots, J$ . Then $p_{J+1} = 1 - \sum_{i=1}^{J} p_{i}$ . Clearly, $\mathbb{P}(N = i) = p_{i}$ for all $i = 1,2,\ldots,J+1$ . Since, any distribution with finite support on $\mathbb{N}$ is a discrete phase-type distribution [Reference Latouche and Ramaswami20], N follows a discrete phase-type distribution. In this case, $N \sim \mathrm{DPH}_{J+1}(\pmb{\mu},\boldsymbol{M})$ such that $\pmb{\mu} = (p_{1},p_{2},\ldots,p_{J+1})$ and

$$\boldsymbol{M} = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots &\ddots & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0\end{array}\right).$$

Shock models with a finite number of shocks were studied in [Reference Montoro-Cazorla and Pérez-Ocón25, Reference Segovia and Labeau34].

6. Numerical illustration

A shock might be seen as a harmful event (e.g. wind gusts for wind turbines, high temperature, etc.). According to the extreme shock model, the failure of a system occurs due to a shock that has a magnitude above a critical threshold. Consider a system which works under the extreme shock model, i.e. a shock can kill the system with probability p and does not affect the system with probability $(1-p)$ . Assume $\mathrm{MAP}(\pmb{\pi},\boldsymbol{D}_{0},\boldsymbol{D}_{1})$ has been found to be useful for modeling the inter-arrival times between subsequent shocks, with parameters

(11) \begin{equation} \pmb{\pi} = (1,0), \quad \boldsymbol{D}_{0} = \left(\begin{array}{c@{\quad}c} -\xi & \xi \\ 0 & -\xi \end{array}\right), \quad \pmb{D}_{1} = \left(\begin{array}{c@{\quad}c} 0 & 0 \\ 0 & \xi \end{array}\right), \qquad \xi > 0.\end{equation}

Thus, with the help of the results derived in Sections 3 and 4, we will be able to compute reliability measures and optimal replacement times for the system. For the extreme shock model, $N \sim \mathrm{DPH}_{1}(\pmb{\mu},\boldsymbol{M})$ where $\pmb{\mu} = 1$ and $\boldsymbol{M} = 1 - p$ . Then, from Theorem 1, $T \sim \mathrm{PH}_{2}(\pmb{\pi}\otimes\pmb{\mu},\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M})$ , where $\pmb{\pi}\otimes\pmb{\mu} = (1,0)$ and

$$\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M} =\left(\begin{array}{c@{\quad}c} -\xi & \xi \\ 0 & -\xi p \\\end{array}\right).$$

Thus, from (6)–(8), we can find expressions for the system’s survival function, the failure rate function, and the mean residual lifetime function. In Figure 1, we plot these reliability indices. Figures 1(a), (c), and (e) show plots for different $p_{i}$ and for the fixed parameter $\xi = 1$ . Figures 1(b), (b), and (f) show plots for different values of $\xi$ . In Figures 1(b) and (d) we fixed parameter $p = 0.2$ , and in Figure 1(f) we fixed parameter $p = 0.12$ . From Figure 1, we can observe that the system’s reliability and mean residual life will decrease if p or $\xi$ increase, while the system’s failure rate increases if p or $\xi$ increase. Clearly, Figures 1(c) and (d) demonstrate the result given in Theorem 3, and $\lim_{t \rightarrow \infty} r_{T}(t) = \xi p$ (the result which we used in Theorem 2).

Figure 1. Plots of the system’s reliability measures.

From (1), $\sigma = \mathbb{E}(T) = ({p+1})/{\xi p}$ . As we can see, the matrix $\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M}$ has distinct eigenvalues for $0 \leq p < 1$ . Therefore this matrix will be diagonalizable, i.e. there exists an invertible matrix $\pmb{Q}$ such that $\boldsymbol{D}_{0}\otimes\boldsymbol{I}+\boldsymbol{D}_{1}\otimes\boldsymbol{M} = \pmb{Q}\,\mathrm{diag}(\!-\!\xi,-\xi p)\pmb{Q}^{-1}$ , where

$$\pmb{Q} = \left(\begin{array}{c@{\quad}c} 1 & 1\\ 0 & 1-p\\\end{array}\right). \qquad \pmb{Q^{-1}} = \frac{1}{(1-p)}\left(\begin{array}{c@{\quad}c} 1-p & -1\\ 0 & 1 \\ \end{array}\right).$$

Note that $(\pmb{\pi}\otimes\pmb{\mu})\boldsymbol{C}_{1}\pmb{R}_{1}\boldsymbol{e} = -p/(1-p) < 0$ for $ 0 < p < 1$ . So, from Theorem 3, $r_{T}(\cdot)$ is strictly increasing for all $0 < p < 1$ . Consequently, from Theorem 2, there exists a finite and unique $x^{*}$ that minimizes the function C(x) in (9) if $\xi p > c_{2}/\sigma (c_{2} - c_{1})$ , equivalently $p > c_{1}/(c_{2} - c_{1})$ , which means an optimal $x^{*}$ exists for all $ c_{1}/(c_{2} - c_{1}) < p < 1$ and $\xi > 0$ . If $0 < p \leq c_{1}/(c_{2} - c_{1})$ then $x^{*} = \infty$ , i.e. a unit should be replaced only at failure. By considering parameters $c_{1} = 1$ , $c_{2} = 10$ , and $\xi = 1$ , we compute the optimal replacement time $x^{*}$ that minimizes (9). According to our derived conditions, an optimal and finite $x^{*}$ exists if $\frac19 < p < 1$ , otherwise $x^{*} = \infty$ . To illustrate this result, we plot the optimal average run cost against x in Figure 2. From Figures 2(a) and (b) we can see that a finite and unique $x^{*}$ exists, since $p = 0.2$ and $p = 0.12$ lie in the interval $\big(\frac19,1\big)$ . On the other hand, from Figures 2(c) and (d) we can see that $x^{*} = \infty$ as $p = \frac19$ and $p = 0.1$ lies in the interval $\big(0,\frac19\big]$ .

Figure 2. Long-run average cost when $c_{1} = 1$ , $c_{2} = 10$ , and $\xi = 1$ .

Now, assume that the system is subject to shocks that occur based on $\mathrm{MAP}(\pmb{\tilde{\pi}},\pmb{\tilde{D}_{0}},\pmb{\tilde{D}_{1}})$ with parameters

(12) \begin{equation} \pmb{\tilde{\pi}} = (1,0), \quad \pmb{\tilde{D}_{0}} = \left(\begin{array}{c@{\quad}c} -\xi & \xi \\ 0 & -\xi \end{array}\right), \quad \pmb{\tilde{D}_{1}} = \left(\begin{array}{c@{\quad}c} 0 & 0 \\ \xi & 0 \end{array}\right), \qquad \xi > 0. \end{equation}

Note that the above process is indeed a phase-type renewal process with inter-arrival times having phase-type representation $(\pmb{\tilde{\pi}},\pmb{\tilde{D}_{0}})$ . The main difference between the Markovian arrival process with parameters given in (11) and the Markovian arrival process with parameters given in (12) is that the former has a dependency between inter-arrival times, and the latter has independent inter-arrival times. However, both processes have the same phase-type generator for inter-arrival times. Figure 3 illustrates the impact of the dependency between the inter-arrival times of the Markovian arrival process on different reliability measures derived in Section 3. In this figure, we compare these reliability measures for the Markovian arrival process in (11) (the process with dependent inter-arrival times) and the phase-type renewal process in (12) (the process with independent inter-arrival times). For Figures 3(a) and (b) we consider the parameters $\xi = 1$ and $p = 0.1$ , while for Figure 33c we consider $\xi = 3$ and $p = 0.2$ . We see that the reliability measures for the dependent and independent inter-arrival times cases are significantly different.

Figure 3. Plots of the system’s reliability measures for phase-type renewal and Markovian arrival shock processes.

For the phase-type renewal process in (12), from Theorem 1, $T \sim \mathrm{PH}_{2}(\tilde{\pmb{\pi}}\otimes\pmb{\mu},\pmb{\tilde{D}_{0}}\otimes\boldsymbol{I}+\pmb{\tilde{D}_{1}}\otimes\boldsymbol{M})$ , where $\tilde{\pmb{\pi}}\otimes\pmb{\mu} = (1,0)$ and

$$\pmb{\tilde{D}_{0}}\otimes\boldsymbol{I}+\pmb{\tilde{D}_{1}}\otimes\boldsymbol{M} =\left(\begin{array}{c@{\quad}c} -\xi & \xi \\ \xi(1 - p) & -\xi \\\end{array}\right).$$

In this case, from (1), $\sigma = 2/\xi p$ and the eigenvalues of the matrix $\pmb{\tilde{D}_{0}}\otimes\boldsymbol{I}+\pmb{\tilde{D}_{1}}\otimes\boldsymbol{M}$ are $-\xi\pm\sqrt{1-p}$ . Therefore, there exists an invertible matrix $\pmb{\tilde{Q}}$ such that $\pmb{\tilde{D}_{0}}\otimes\boldsymbol{I}+\pmb{\tilde{D}_{1}}\otimes\boldsymbol{M} = \pmb{\tilde{Q}}\,\mathrm{diag}(-\xi - \xi\sqrt{1-p},-\xi + \xi\sqrt{1-p})\,\pmb{\tilde{Q}}^{-1}$ , where

$$\pmb{\tilde{Q}} = \left(\begin{array}{c@{\quad}c} 1 & 1\\ -\sqrt{1-p} & \sqrt{1-p}\\\end{array}\right), \qquad \pmb{\tilde{Q}}^{-1} = \frac{1}{2\sqrt{1-p}}\left(\begin{array} {c@{\quad}c} \sqrt{1-p} & -1\\ \sqrt{1-p} & 1 \\\end{array}\right). $$

Clearly, $(\tilde{\pmb{\pi}}\otimes\pmb{\mu})\boldsymbol{C}_{1}\pmb{R}_{1}\boldsymbol{e} = (\sqrt{1-p} - 1)/2(\sqrt{1-p}) < 0$ for $ 0 < p < 1$ . So, from Theorem 3, $r_{T}(\cdot)$ is strictly increasing for all $0 < p < 1$ . Consequently, from Theorem 2, there exists a finite and unique $x^{*}$ that minimizes the function C(x) in (9) if $\xi - \xi\sqrt{1-p} > c_{2}/\sigma (c_{2} - c_{1})$ , equivalently $p > 4 c_{1} (c_{2} - c_{1})/c_{2}^{2}$ , which means an optimal $x^{*}$ exists for all $4c_{1}(c_{2} - c_{1})/c_{2}^{2} < p < 1$ and $\xi > 0$ . If $0 < p \leq 4c_{1}(c_{2} - c_{1})/c_{2}^{2}$ then $x^{*} = \infty$ . Let $c_{1} = 1$ and $c_{2} = 10$ . In this case, if $0.36 < p < 1$ then a finite and unique $x^{*}$ exists, and the corresponding per unit mean cost is given by (9). Observe that dependency between the inter-arrival times can change decisions significantly. For example, for $p = 0.2$ , in the case of the Markovian arrival process in (11), a finite $x^{*}$ exists, i.e. if the system does not fail by time $x^{*}$ , then we should replace our system with a new one at time $x^{*}$ . However, in the renewal process in (12), $x^{*} = \infty$ , i.e. we should replace the system only at failure. Even if an optimal $x^{*}$ exists for both cases, the values of $x^{*}$ and $C(x^{*})$ differ significantly in both cases. Therefore, consideration of dependency between inter-arrival times can improve our decision-making.

Let $x_{1}^{*}$ (respectively $C_{1}(x_{1}^{*})$ ) and $x_{2}^{*}$ (respectively $C_{2}(x_{2}^{*})$ ) denote the optimal values of the replacement time (respectively the optimal value of the mean cost function) for the Markovian arrival process with parameters given in (11) and (12), respectively. Table 1 shows these optimal values for different values of p and $\xi$ , which shows that the values of $x_{1}^{*}$ (when inter-arrival times are dependent) and the values of $x_{2}^{*}$ (when inter-arrival times are independent) are significantly different, and correspondingly for the optimal mean cost. Further, we observe that an increment in the value of p or $\xi$ decreases the optimal replacement time and increases the value of the optimal mean cost (irrespective of the dependency between inter-arrival times).

Table 1. Optimal values of x and C(x) for different values of p and $\xi$ .