A general class of shock models with dependent inter-arrival times

Abstract

We introduce and study a general class of shock models with dependent inter-arrival times of shocks that occur according to the homogeneous Poisson generalized gamma process. A lifetime of a system affected by a shock process from this class is represented by the convolution of inter-arrival times of shocks. This class contains many popular shock models, namely the extreme shock model, the generalized extreme shock model, the run shock model, the generalized run shock model, specific mixed shock models, etc. For systems operating under shocks, we derive and discuss the main reliability characteristics (namely the survival function, the failure rate function, the mean residual lifetime function and the mean lifetime) and study relevant stochastic comparisons. Finally, we provide some numerical examples and illustrate our findings by the application that considers an optimal mission duration policy.

1 Introduction

Shock models are widely used for stochastic description of systems operating in random environments. Various shock models were introduced and studied in the literature. One can classify most of the shock models into four major groups, namely the extreme shock model (see, e.g., Gut and Hüsler 1999, 2005; Shanthikumar and Sumita 1983, 1984; Cha and Finkelstein (2016)), the cumulative shock model (see, e.g., A-Hameed, M.S. and Proschan, F. (1973); Esary et al. (1973); Gut (1990); Gong et al. (2020a)), the run shock model (see, e.g., Mallor and Omey (2001); Gong et al. (2018)) and the $\delta$-shock model (see, e.g., Li and Kong (2007); Goyal et al. (2022b, 2022c)). For convenience, we provide definitions of these models in Subsect. 2.3. Besides, there are various mixed shock models that are combinations of two or more shock models from these major groups (see Cha and Finkelstein (2009); Eryilmaz and Tekin (2019); Wang and Zhang (2005); Mallor et al. (2006); Eryilmaz (2012), and the references therein).

Most of the studies devoted to shock modeling consider renewal processes or Poisson processes (homogeneous Poisson process (HPP) and non-homogeneous Poisson process (NHPP)) as the corresponding shock point processes affecting engineering systems. Note that the inter-arrival times in the renewal process are i.i.d., which is very restrictive in many real-life scenarios, whereas the Poisson process has independent increments, which is also a restrictive assumption in many applications. For example, if there are larger number of shocks in the past, we may expect the same in the future (positive dependence). Thus, in many real-life scenarios, it is more appropriate to consider more general shock process.

Recently, a new point process called the Poisson generalized gamma process (PGGP) has been introduced by Cha and Mercier (2021). It has a complex mathematical structure, which makes its implementation in applications limited. However, its special case, namely the homogeneous Poisson generalized gamma process (HPGGP), not only is mathematically tractable but also has many nice statistical properties. The HPGGP possesses the dependent increments property with a general dependency structure and hence, it may be applied to a wider class of problems. Further, the inter-arrival times of this process are also dependent. Moreover, it contains many well-known processes (namely the HPP, the Pólya process, etc.) as the particular cases. Thus, the HPGGP can be viewed as a more general and flexible process that does not have limitations of the Poisson and the renewal processes.

It should be noted that the Pólya process, being a special case of the HPGGP, is also widely used in shock modeling. The inter-arrival times of this process follow the Pareto distribution that is well known in modeling events with heavy-tailed behavior. For example, an earthquake or a storm can be considered as a shock for some critical systems, e.g., a bridge. If the magnitude of an earthquake is sufficiently large, then the bridge can be severely damaged or completely destroyed (see, e.g., Last and Szekli (1998); Eryilmaz (2017b); Cha and Finkelstein (2016) based on the Pólya shock process). Some of the other specific settings were considered in our recent works (Goyal et al. (2022a, 2022b, 2022c)).

On the other hand, there exists many popular shock models, which have a great potential in different applications, namely the run shock model, the generalized run shock model, the generalized extreme shock model, and some special mixed shock models. To the best of our knowledge, they have not yet been studied with sufficient generality taking into account the dependent inter-arrival times. Thus, in this paper, our goal is to develop and study a unified class of shock models with dependent inter-arrival times of shocks that occur according to the HPGGP. Moreover, we assume that the distribution of the fatal shock that causes the system’s failure (i.e., the number of shocks until the system failure) is a discrete phase-type (DPH) distribution, which describes a rather general model of failure. The survival function and other statistical measures (e.g., mean, variance, etc.) of the DPH can be written in matrix forms that are convenient in computations by using the relevant mathematical packages. Consequently, the results to be developed for the proposed class of shock models are also mathematically tractable.

The contribution and novelty of our paper can be summarized as follows:

(a):

We propose the unified class of shock models with the specified dependence structure and a general discrete failure model. This class contains many popular shock models, namely the run shock model, the generalized run shock model, the generalized extreme shock model, some special mixed shock models, etc.;

(b):

We assume that shocks occur in accordance with the homogeneous Poisson generalized gamma process (HPGGP), which is a rather general process with dependent increments and inter-arrival times, whereas the failure model (i.e., the distribution of the number of shocks until the system failure) is defined by the discrete phase-type (DPH) distribution;

(c):

For the proposed class of shock models, we study some important reliability indices, namely the survival function, the failure rate function, the mean residual life function and the mean lifetime. Moreover, the derived results are in matrix forms that can easily be calculated using the relevant mathematical packages.

The rest of the paper is organized as follows. In Sect. 2, we provide some preliminaries and then define the class of shock models to be considered. In Sect. 3, we derive some important reliability indices and study some stochastic comparisons for the proposed class of shock models. In Sect. 4, we give some numerical examples to illustrate the developed results. In Sect. 5, we illustrate our findings by considering the relevant optimal mission duration policy. Lastly, the concluding remarks are given in Sect. 6.

To enhance the readability of the paper, all proofs of theorems, lemmas and corollaries, wherever given, are deferred to the Appendix, which forms a supplementary material to our paper.

2 Preliminaries and model formulation

For any random variable U, we denote the cumulative distribution function (cdf) by $F_U(\cdot )$, the survival/reliability function by ${\bar{F}}_U(\cdot )$, the probability density function (pdf) by $f_{U}(\cdot )$ and the failure rate function by $r_{U}(\cdot )$; here ${\bar{F}}_U(\cdot )\equiv 1-F_U(\cdot )$ and $r_{U}(\cdot ) \equiv f_{U}(\cdot )/\bar{F}_{U}(\cdot )$. Further, we denote the set of natural numbers by ${\mathbb {N}}$, the set of real numbers by ${\mathbb {R}}$ and the set of complex numbers by ${\mathbb {C}}$. For any $z \in {\mathbb {C}}$, Re(z) and Im(z) denote the real part and the imaginary parts of z, respectively. By writing a matrix ${\varvec{A}} = diag(x_{1},x_{2},\dots ,x_{d})$, we mean that ${\varvec{A}}$ is diagonal matrix of size d with ith diagonal entry equal to $x_{i}$, $i = 1,2,\dots ,d$.

2.1 Shock processes

Let $\{N(t),\;t\ge 0\}$ be an orderly point process representing the number of shocks arrived by time t. In the literature, different point processes of shocks were developed to model different real-life scenarios (see Cha and Finkelstein (2018), Teugels and Vynckier (1996), and the references therein). In what follows, we give the formal definitions of some point processes that will be used in this paper. But first, we recall the definition of the generalized gamma distribution (see Agarwal and Kalla (1996)).

Definition 2.1

A random variable Q is said to have the generalized gamma distribution (GGD) with the set of parameters $\{\nu ,k,\alpha ,l\}$, $\nu \ge 0,\; k,\alpha ,l > 0$, denoted by $Q \sim GG(\nu ,k,\alpha ,l)$, if its pdf is given by

$$\begin{aligned} f_{Q}(q) = \frac{\alpha ^{k-\nu }}{\Gamma _{\nu }(k,\alpha l)} \frac{q^{k-1} \exp \{ - \alpha q \}}{(q + l)^{\nu }}, \quad q > 0, \end{aligned}$$

(1)

where

$$\begin{aligned} \Gamma _{\nu }(k,\alpha l) = \int _{0}^{\infty } \frac{y^{k-1} \exp \{ -y \}}{( y + \alpha l )^{\nu }} dy = \int _{0}^{\infty } \frac{ \alpha ^{k - \nu } y^{k-1} \exp \{ - \alpha y \}}{( y + l )^{\nu }} dy. \end{aligned}$$

(2)

Definition 2.2

A mixed Poisson process $\{ N(t), t \ge 0 \}$ is said to be the Pólya process with the set of parameters $\{\beta ,b\},\; \; \beta>0,\;b>0$, if

$$\begin{aligned} \begin{aligned} P(N(t) = n)&= \frac{\Gamma (\beta + n)}{\Gamma (\beta ) n!} \left( \frac{t}{t+b} \right) ^{n} \left( \frac{b}{t+b} \right) ^{\beta } \\&= \int _{0}^{\infty } \exp \{ - \chi t\} \frac{(\chi t )^{n} }{ n !} dH(\chi ),\quad n=0,1,2,\dots , \end{aligned} \end{aligned}$$

where H, named as the structure distribution, is the gamma distribution with the density given by

$$\begin{aligned} dH(\chi ) = \frac{b^{\beta }}{\Gamma (\beta )} \chi ^{\beta -1} \exp \{ - b \chi \} d \chi . \end{aligned}$$

Remark 2.1

Let $X_{1}, X_{2}, \dots$ be a sequence of inter-arrival times of the Pólya process with set of parameters $\{\beta , b\}$. Then, $X_{i}$ follows the Pareto distribution with the cdf given by

$$\begin{aligned} F_{X_{i}}(t)= 1 - \left( \frac{b}{b + t}\right) ^{\beta },\;\;t\ge 0,\;\;i=1,2,3,\dots \end{aligned}$$

Definition 2.3

A counting process $\{ N(t), t \ge 0 \}$ is said to be the Poisson generalized gamma process (PGGP) with the set of parameters $\{\lambda (t),\nu ,k,\alpha ,l\}, \; \lambda (t) > 0 \text { for all } t \ge 0, \; \nu \ge ~0,$ $k,\alpha ,l >~0$, if

(a):

$\{ N(t), t \ge 0 |Q = q\} \sim NHPP(q \lambda (t))$;

(b):

$Q \sim GG(\nu ,k,\alpha ,l)$.

If $\lambda (t) = \lambda$ $(>0)$, then it is called the homogeneous Poisson generalized gamma process (HPGGP) with the set of parameters $\{\lambda , \nu , k,\alpha , l\}$. $\square$

Note that the PGGP is a mixed Poisson process with the generalized gamma as a mixing distribution. Moreover, it can be understood as a Poisson process having the stochastic intensity function which is a product of the deterministic intensity function and a random variable.

As we discussed in the Introduction, the HPGGP possesses the dependent increments property. The dependence structure in the increments of the HPGGP is given by the positive upper orthant dependent increments property (see Cha and Mercier (2021)), i.e., for any arbitrary integer $m \ge 2$ and $0< t_{1}< t_{2}< \dots < t_{m}$,

$$\begin{aligned} P(N(t_{i} + \Delta t_{i}) - N(t_{i})> n_{i}, i = 1,2,\dots , m) \ge \prod _{i=1}^{m} P(N(t_{i} + \Delta t_{i}) - N(t_{i}) > n_{i}), \end{aligned}$$

for all $n_{i}$, $i = 1,2,\dots ,m$, where $t_{i} + \Delta t_{i} \le t_{i+1}$, $i = 1, 2, \dots , m-1$.

Remark 2.2

The following statements are true

(a):

The HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$, where $\nu = 0,\alpha =k$ and $k\rightarrow \infty$, is the HPP with the intensity $\lambda$, regardless of l;

(b):

The HPGGP with the set of parameters $\{\lambda ,\nu ,k,\alpha ,l\}$, where $\lambda = 1/b$ $(>0)$, $\nu = 0, k = \beta , \alpha = 1$, is the Pólya process with the set of parameters $\{\beta ,b\}$, regardless of l.

2.2 Discrete PH distribution

Below, we give the definition of a discrete PH distribution (see, e.g., Eryilmaz (2017a), Bozbulut and Eryilmaz (2020)).

Definition 2.4

A discrete phase-type distribution can be viewed as the distribution of the time to absorption in an absorbing Markov chain. A random variable N is said to have the discrete PH distribution, denoted by $N \sim DPH({\varvec{a}}, {\varvec{A}})$, if

$$\begin{aligned} P(N = n) = {\varvec{a}} {\varvec{A}}^{n-1} {\varvec{u}},\quad n \in {\mathbb {N}}, \end{aligned}$$

where ${\varvec{A}}$ is a square matrix of size m that includes the transition probabilities among the m transient states, and ${\varvec{a}} = (a_{1},a_{2},\dots ,a_{m})$ is a vector of size m such that $\sum _{i=1}^{m} a_{i} = 1$. Moreover, ${\varvec{u}} = ({\varvec{I}} - {\varvec{A}}){\varvec{e}}$ is a vector which includes the transition probabilities from transient states to the absorbing state, ${\varvec{a}} {\varvec{e}} = 1$, and ${\varvec{I}}$ is the identity matrix, and ${\varvec{e}}$ is the column vector of size m with all elements being one. Further, the matrix ${\varvec{A}}$ must satisfy the condition that ${\varvec{I}} - {\varvec{A}}$ is non-singular. $\Box$

    Note that the survival function of N is given by $\bar{F}_{N}(n) = {\varvec{a}} {\varvec{A}}^{n} {\varvec{e}}$ and $E(N) = {\varvec{a}} ({\varvec{I}} - {\varvec{A}})^{-1} {\varvec{e}}$.

2.3 Model formulation and specific cases

Let L be a random variable defining a lifetime of a system incepted into operation at time $t= 0$ and subject to external random shocks. We assume that shocks is the only cause of its failure. Let $X_{i}$ be a random variable representing the time between the $(i-1)$th and the ith shocks, $i \ge 1$. Further, let $Y_{i}$ be a random variable defining the magnitude (or, the damage caused by) by the ith shock, $i = 1,2,3,\dots$. Furthermore, let M be a random variable representing the fatal shock, i.e., the number of shocks arrived till the failure of the system. We define now a general class of shock models for which the lifetime of a system can be represented as a random sum of the inter-arrival times of shocks, i.e., $L = \sum _{i=1}^{M} X_{i}$. We study this class of shock models under the following assumptions.

Assumptions:

• Shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$;

• $M \sim DPH({\varvec{a}},{\varvec{A}})$;

• $\{Y_{i}: i \in {\mathbb {N}}\}$ is a sequence of independent and identically distributed random variables. Moreover, $\{Y_{i}: i \in {\mathbb {N}}\}$ and $\{X_{i}: i \in {\mathbb {N}}\}$ are independent;

• M and $\{X_{i}: i \in {\mathbb {N}}\}$ are independent.

Below, we list some well-known shock models that belong to this class.

• Extreme shock model (Gut and Hüsler 1999, 2005): In this model, a system fails at the ith shock if its magnitude exceeds the predetermined threshold $\gamma$, i.e., $Y_{i} > \gamma$. Consequently,

  $$\begin{aligned} P(M = m)= & {} P(Y_{1} \le \gamma , Y_{2} \le \gamma , \dots , Y_{m-1} \le \gamma , Y_{m} > \gamma ) = p(1-p)^{m-1}, \end{aligned}$$

  where $p = P(Y_{1} > \gamma )$. Clearly, $M \sim DPH(1,1 - p)$. Further, note that M and $\{X_{i}: i \in {\mathbb {N}}\}$ are independent. Thus, the extreme shock model is a member of the defined class of shock models.

• Generalized extreme shock models I and II (Bozbulut and Eryilmaz (2020)): In these models, it is assumed that shocks may arrive from more than one possible sources (say, n sources). Let $\theta _{i}$ be the probability that shocks come from source i, and let $p_i$ be the probability that the magnitude of a shock from source i exceeds the critical level $\gamma$; $i = 1,2,\dots ,n$. In both models, a system fails upon occurrence of a shock of size, at least $\gamma$. Thus, shocks that come from source i are harmless with probability $1 - p_{i}$; $i = 1,2,\dots ,n$. In Model I, it is assumed that shocks may arrive from different sources over time, i.e., two consecutive shocks may come from two different sources. Note that, in this model, $M \sim DPH(1, 1-p)$, where $p = \sum _{i = 1}^{n} \theta _{i} p_{i}$ and $\sum _{i=1}^{n} \theta _{i} = 1$ (see Bozbulut and Eryilmaz (2020)). In Model II, it is assumed that shocks initially come from a specified source, say i, and after the ith source completes its task, shocks continue to arrive from another source. In this case, $M \sim DPH( {\varvec{a}_{gem}}, {\varvec{A}_{gem}})$, where ${\varvec{a}_{gem}} = (\theta _{1}, \theta _{2},\dots ,\theta _{n})$, ${\varvec{A}_{gem}} = diag(1-p_{1}, 1-p_{2}, \dots , 1-p_{n})$ and $\sum _{i=1}^{n} {\theta _{i}} = 1$. If $n = 1$, then both models reduce to the classical extreme shock model. Clearly, in both models, M follows the discrete PH distribution, and M and $\{X_{i}: i \in {\mathbb {N}}\}$ are independent. Therefore, these two models belong to the defined class of shock models.

• Run shock model (Mallor and Omey (2001)): In this model, a system fails when the magnitudes of r consecutive shocks exceed the predetermined threshold value $\gamma$. Here, the random variable M follows the geometric distribution of order r with the success probability p, where p is the probability that the magnitude of a single shock exceeds the threshold value $\gamma$. Later, Tank and Eryilmaz (2015) showed that the distribution of M can equivalently be represented as a discrete PH distribution i.e., $M \sim DPH( {\varvec{a_{rm}}}, {\varvec{A_{rm}}})$, where ${\varvec{a_{rm}}} = (1,0,\dots ,0)_{1 \times r}$ and

  $$\begin{aligned} {\varvec{A_{rm}}} = \begin{pmatrix} 1 - p &{}\quad p &{}\quad 0 &{}\quad \cdots &{}\quad 0 \\ 1 - p &{}\quad 0 &{}\quad p &{}\quad \cdots &{}\quad 0 \\ \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ 1-p &{}\quad 0 &{}\quad 0 &{}\quad \cdots &{}\quad p\\ 1 - p&{}\quad 0 &{}\quad 0 &{}\quad \cdots &{}\quad 0 \end{pmatrix}_{r \times r}. \end{aligned}$$

  Consequently, the run shock model belongs to the defined class of shock models.

• Mixed shock model (Eryilmaz and Tekin (2019)): In this model, the classical run and extreme shock models are combined. Let $\gamma _{1}$ and $\gamma _{2}$ be two fixed threshold values such that $\gamma _{1} < \gamma _{2}$. Here, the failure of a system takes place when either r consecutive shocks of size, at least $\gamma _{1}$ or a single shock of size, at least $\gamma _{2}$ occur(s). In this case, $M \sim DPH({\varvec{a_{mm}}}, {\varvec{A_{mm}}})$, where ${\varvec{a_{mm}}} = (1,0,\dots ,0)_{1\times r}$ and

  $$\begin{aligned} {\varvec{A_{mm}}} = \begin{pmatrix} 1 - p_{1} - p_{2} &{} p_{1} &{}0 &{} \cdots &{} 0 \\ 1 - p_{1} - p_{2} &{} 0 &{} p_{1} &{} \cdots &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots \\ 1 - p_{1}- p_{2} &{}0 &{} 0 &{} \cdots &{} p_{1} \\ 1 - p_{1}- p_{2} &{}0 &{} 0 &{} \cdots &{} 0 \end{pmatrix}_{r \times r}, \end{aligned}$$

  with $p_{1} = P(\gamma _{1} \le Y_{i} < \gamma _{2})$ and $p_{2} = P(Y_{i} \ge \gamma _{2})$ for $i = 1,2,\dots$. Consequently, this model also belongs to the defined class of shock models.

• Generalized run shock model I (Gong et al. (2018)): Consider two critical levels $\gamma _1$ and $\gamma _2$ such that $\gamma _{1} < \gamma _{2}$. For given two positive integers $r_{1}$ and $r_{2}$, the system is assumed to fail if, at least $r_{1}$ consecutive shocks with magnitude above $\gamma _{1}$ or $r_{2}$ consecutive shocks with magnitude above $\gamma _{2}$ occur. It is obvious that $r_{1} > r_{2}$. Let $p_{1} = P(Y_{i} \le \gamma _{1})$, $p_{2} = P(\gamma _{1} < Y_{i} \le \gamma _{2})$ and $p_{3} = P(Y_{i} > \gamma _{2})$, $i = 1,2,3,\dots$. Then, $M \sim DPH({\varvec{a_{grm1}}}, {\varvec{A_{grm1}}})$, where ${\varvec{a_{grm1}}} = (1, 0, \dots , 0)$ and

  $$\begin{aligned} {\varvec{A_{grm1}}}= & {} \begin{pmatrix} {\varvec{E_{r_{1}}}} + {\varvec{S_{r_{1}}}} &{} \quad {\varvec{T_{r_{1}}}} &{} \quad \\ {\varvec{E_{r_{1}-1}}} &{}\quad {\varvec{S_{r_{1}-1}}} &{}\quad {\varvec{T_{r_{1}-1}}} \\ {\varvec{E_{r_{1}-2}}} &{} &{}\quad {\varvec{S_{r_{1}-2}}} &{} \quad \ddots \\ \vdots &{} &{} &{}\quad \ddots &{}\quad {\varvec{T_{r_{1}-r_{2} + 2}}}\\ {\varvec{E_{r_{1} - r_{2} +1}}} &{} &{} &{} &{} {\varvec{S_{r_{1}-r_{2} + 1}}} \end{pmatrix}, \\ {\varvec{E_{j}}}= & {} \begin{pmatrix} p_{1} &{}\quad 0 &{}\quad \dots &{}\quad 0\\ p_{1} &{}\quad 0 &{}\quad \dots &{}\quad 0\\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ p_{1} &{}\quad 0 &{} \quad \dots &{}\quad 0 \end{pmatrix}_{j \times r_{1}}, \\ {\varvec{S_{j}}}= & {} \begin{pmatrix} 0 &{}\quad p_{2} &{}\quad \dots &{}\quad 0 \\ 0 &{}\quad \vdots &{}\quad \ddots &{}\quad 0 \\ \vdots &{}\quad 0 &{}\quad \dots &{}\quad p_{2} \\ 0 &{} \quad 0 &{}\quad \dots &{} \quad 0 \end{pmatrix}_{j \times j}, \quad {\varvec{T_{j}}} = \begin{pmatrix} p_{3} &{}\quad \dots &{}\quad 0\\ \vdots &{}\quad \ddots &{}\quad 0\\ 0 &{}\quad \dots &{}\quad p_{3} \\ 0 &{}\quad \dots &{}\quad 0 \end{pmatrix} _{j \times j-1}. \end{aligned}$$

  Clearly, this model belongs to the defined class.

• Generalized run shock model II (Yalcin et al. (2018)): In this model, a system fails due to the occurrences of n non-overlapping runs of r consecutive critical shocks. This model can be useful for describing the performance of the n cold standby identical systems because, in this case, a system needs to be replaced with a new one after the occurrence of r consecutive critical shocks. For this model, $M \sim DPH( {\varvec{a_{grm2}}}, {\varvec{A^{(n)}_{grm2})}}$, where

  $$\begin{aligned} {\varvec{a_{grm2}}} = (1,0,\dots ,0)_{1\times nr}, \quad {\varvec{A^{(n)}_{grm2}}} = \begin{pmatrix} {\varvec{A^{(n-1)}_{grm2}}} &{} {\varvec{S}} \\ 0 &{} {\varvec{T}}\end{pmatrix}, \quad n \ge 3, \end{aligned}$$

  and ${\varvec{A^{(n)}_{grm2}}}$ is an $nr \times nr$ matrix with

  $$\begin{aligned} {\varvec{A^{(2)}_{grm2}}} = \begin{pmatrix} 1-p &{} \quad p &{} \quad 0 &{}\quad \dots &{}\quad 0 &{}\quad 0 &{}\quad 0 &{} \quad \dots &{}\quad 0 \\ 1-p &{} 0 &{} p &{} \dots &{} 0 &{} 0 &{} 0 &{} \dots &{} 0 \\ \vdots &{}\quad \vdots &{} \quad \vdots &{}\quad \vdots &{} \quad \vdots &{}\quad \vdots &{} \quad \vdots &{}\quad \vdots &{}\quad \vdots \\ 1-p &{} \quad 0 &{}\quad 0 &{}\quad \dots &{}\quad 0 &{}\quad p &{} \quad 0 &{}\quad \dots &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{} \quad \dots &{}\quad 0 &{}\quad 1-p &{}\quad p &{}\quad \dots &{} \quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots &{}\quad 0 &{}\quad 1-p &{}\quad 0 &{}\quad \dots &{}\quad 0 \\ \vdots &{} \quad \vdots &{}\quad \vdots &{}\quad \vdots &{} \quad \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \vdots \\ 0 &{} \quad 0 &{}\quad 0 &{}\quad \dots &{}\quad 0 &{}\quad 1-p &{}\quad 0 &{} \quad \dots &{}\quad 0 \\ \end{pmatrix}_{2r \times 2r}. \end{aligned}$$

  Further, the matrices ${\varvec{S}}$ and ${\varvec{T}}$ are given by

  $$\begin{aligned} {\varvec{S}} = \begin{pmatrix} 0 &{}\quad 0 &{}\quad 0 &{} \quad \dots &{}\quad 0 \\ \vdots &{}\quad \vdots &{}\quad \vdots &{} \quad \dots &{}\quad \vdots \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad \vdots &{} \quad 0 \\ p &{}\quad 0 &{}\quad 0 &{} \quad \dots &{}\quad 0 \end{pmatrix}_{(n-1)r \times r}, \quad {\varvec{T}} = \begin{pmatrix} 1-p &{}\quad p &{}\quad 0 &{}\quad \dots &{}\quad 0 \\ 1-p &{}\quad 0 &{}\quad p &{} \quad \dots &{} \quad 0 \\ \vdots &{} \quad \vdots &{}\quad \vdots &{}\quad \dots &{}\quad \vdots \\ 1-p &{}\quad 0 &{}\quad 0 &{}\quad \dots &{}\quad 0 \end{pmatrix}_{r \times r}. \end{aligned}$$

  Thus, this model also belongs to our class.

The suggested in our paper class of shock models is general, but obviously, it does not contain all existing models. The following examples illustrate this reasoning:

• Cumulative shock model (A-Hameed, M.S. and Proschan, F. (1973)): In this model, a system fails at ith shock if the cumulative damage till the occurrence of ith shock exceeds the predetermined threshold value $\gamma$ and the cumulative damage till the occurrence of the $(i-1)$th shock is less than $\gamma$, i.e., $\sum _{j=1}^{i} Y_{j} > \gamma$ and $\sum _{j=1}^{i-1} Y_{j} \le \gamma$. Then,

  $$\begin{aligned} P(M = m) = P\left( \sum _{j=1}^{m-1} Y_{j} \le \gamma , \sum _{j=1}^{m} Y_{j} > \gamma \right) = F^{(m-1)}(\gamma ) - F^{(m)}(\gamma ), \end{aligned}$$

  where $F^{(m)}(\cdot )$ is the cdf of $\sum _{j=1}^{m}Y_{j}$. In general, M does not follow the discrete PH distribution and hence, the cumulative shock model is not a member of the defined class of shock models. However, if the distribution of M follows the discrete PH distribution, for a particular choice of $F^{(m)}(\cdot )$, then the cumulative shock model will also be a member of the defined class of shock models.

• $\delta$-shock model (Li and Kong (2007)): In this model, a system fails when the time between two consecutive shocks is less than the predetermined threshold $\delta$. Here, $P(M = m) = P(X_{1}> \delta , X_{2}> \delta , \dots , X_{m-1} > \delta , X_{m} \le \delta ).$ Clearly, M and $\{X_{i}: i \in {\mathbb {N}}\}$ are not independent. Thus, the $\delta$-shock model does not belong to our class.

3 Survival of systems under shocks

In this section, we discuss some reliability indices (namely the survival function, the failure rate function, the mean residual lifetime function and the mean lifetime) for the defined class of shock models. We also discuss some relevant stochastic comparisons.

3.1 Survival function

We start with three lemmas that will be used in proving the main results of this subsection. The proof of the first lemma follows from the fact that the inter-arrival times of the HPP are i.i.d. and follow the exponential distribution (see also equation (6) of Eryilmaz (2017a)).

Lemma 3.1

Let $M \sim DPH({\varvec{a}}, {\varvec{A}})$. Assume that shocks occur according to the HPP with intensity $\theta > 0$. Then, the corresponding survival function for this model is given as $\bar{F}_{L}(t) = {\varvec{a}} \exp \{ - \theta t ({\varvec{I}} - {\varvec{A}}) \} {\varvec{e}}.$ $\Box$

    Before stating the next lemma, we give the formal definition of the Fourier transform of a function (see Osgood (2019), p. 105).

Definition 3.1

Let $g: {\mathbb {R}} \rightarrow {\mathbb {R}}$. Then, the Fourier transform of g, denoted by ${\mathcal {F}}[g](s)$, is given by

$$\begin{aligned} {\mathcal {F}}[g](s) = \int _{-\infty }^{\infty } \exp \{- 2 \pi i s t\} g(t) dt, \text { for }s \in {\mathbb {R}}, \end{aligned}$$

provided that the integral exists.

Lemma 3.2

(Osgood (2019), pp. 114,123 and 135) The following results hold true.

• (i) Let $g(t) = \exp \{- c t\} I_{0}(t)$, where

  $$\begin{aligned} I_{0}(t) = {\left\{ \begin{array}{ll} 0 &{} t\le 0 \\ 1 &{} t > 0 \end{array}\right. } \end{aligned}$$

  and c is a positive constant. Then, ${\mathcal {F}}[g](s) = {1}/{(2 \pi i s + c)};$

• (ii) ${\mathcal {F}}[c g](s) = c {\mathcal {F}}[g](s)$, where c is a constant number (real or complex);

• (iii) $\frac{d^{k}}{ds^{k}} {\mathcal {F}}[ g](s) = {\mathcal {F}}[(- 2 \pi i t)^{k} g(t)](s).$

Lemma 3.3

Let $\epsilon \in {\mathbb {C}}$ such that $Re(\epsilon ) > 0$. Then, for any $n \in {\mathbb {N}}$,

$$\begin{aligned} \int _{0}^{\infty }\exp \{-(\epsilon x +1) q \} q^{n-1} dq = \Gamma (n) \left( \frac{1}{1+\epsilon x} \right) ^{n}, \quad x \ge 0. \end{aligned}$$

     The proof of this lemma and of all other results, for convenience, is deferred to the Appendix (which forms supplementary materials in our paper). Although some of our results look cumbersome, their structure is probabilistically rather simple. Besides, they are given in the matrix form that is useful for computations in applications (see later). In the following theorem, the survival function of a system under the considered shock process is derived.

     Note that there are no assumptions on the structure of matrix ${\varvec{A}}$ (see Lemma 3.1) for the shock processes with independent inter-arrival times (e.g., the HPP). However, when the shock processes with dependent inter-arrival times (e.g., PGGP/Polya process) are considered, some sufficient conditions on the eigenvalues of the matrix ${\varvec{A}}$ (e.g., every eigenvalue of the matrix ${\varvec{I}} - {\varvec{A}}$ is real and positive, real part of every eigenvalue of the matrix ${\varvec{I}} - {\varvec{A}}$ is positive) are implemented. These conditions are feasible and easy to verify. Moreover, they hold for aforementioned specific shock models (see examples in Section 4).

Theorem 3.1

Let $M \sim DPH({\varvec{a}},{\varvec{A}})$, where ${\varvec{A}}$ is a square matrix of size d. Further, let $\epsilon _{1},\epsilon _{2},\dots ,\epsilon _{r}$ be eigenvalues of the matrix ${\varvec{I}}-{\varvec{A}}$ with the algebraic multiplicity $d_{1},d_2,\dots ,d_r$, respectively, and let ${\varvec{P}}$ be a non-singular matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag({\varvec{J_{1}}}, {\varvec{J_{2}}}, \dots , {\varvec{J_{r}}}) {\varvec{P^{-1}}}$, where ${\varvec{J_{i}}}$ is the $\epsilon _{i}$-Jordan block with dimension $d_{i}$, $i = 1,2,\dots ,r$, and $\sum _{i=1}^{r}d_{i} = d$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. If every eigenvalue of the matrix ${\varvec{I}} - {\varvec{A}}$ is real and positive, then the corresponding survival function for a system under the considered shock process is

  $$\begin{aligned} \bar{F}_{L}(t)= & {} \left( \frac{\alpha ^{k-\nu }}{\Gamma _{\nu }(k,\alpha l)} \right) {\varvec{a}} {\varvec{P}} diag({\varvec{D_{1}(t)}}, {\varvec{D_{2}(t)}}, \dots , {\varvec{D_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}, \end{aligned}$$

  (3)

  where

  $$\begin{aligned} {\varvec{D_{i}(t)}}= & {} \begin{pmatrix} \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k-\nu }} &{} \frac{ - \lambda t\Gamma _{\nu }(k + 1, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k + 1-\nu }} &{} \dots &{} \frac{(-\lambda t)^{d_{i}-1} \Gamma _{\nu }(k + d_{i} -1, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-1)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i} -1-\nu }} \\ 0 &{} \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k-\nu }} &{} \dots &{} \frac{(-\lambda t)^{d_{i}-2} \Gamma _{\nu }(k + d_{i} -2, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-2)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i} -2-\nu }}\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k-\nu }} \end{pmatrix}_{d_{i}\times d_{i}}, \end{aligned}$$

  for $t \ge 0$, $i = 1,2,\dots ,r$.

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,r$, then

  $$\begin{aligned} \bar{F}_{L}(t) = \frac{1}{\Gamma (\beta )} {\varvec{P}} diag({\varvec{{\widetilde{D}}_{1}(t)}}, {\varvec{{\widetilde{D}}_{2}(t)}}, \dots , {\varvec{{\widetilde{D}}_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}, \end{aligned}$$

  where

  $$\begin{aligned} {\varvec{{\widetilde{D}}_{i}(t)}} = \begin{pmatrix} \Gamma (\beta ) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta } &{} - \Gamma (\beta + 1) \left( \frac{t}{t \epsilon _{i} + b}\right) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta }&{} \dots &{} \frac{\Gamma (\beta + d_{i}-1)}{(d_{i}-1)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-1} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta } \\ 0 &{} \Gamma (\beta ) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta } &{} \dots &{} \frac{\Gamma (\beta + d_{i}-2)}{(d_{i}-2)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-2} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta }\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \Gamma (\beta ) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta } \end{pmatrix}_{d_{i}\times d_{i}}, \end{aligned}$$

  for $t \ge 0$, $i = 1,2,\dots ,r$. Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,r$, then the above result holds for all $\beta \in {\mathbb {R}}$.$\Box$

The following corollary immediately follows from Theorem 3.1, Remark 2.2(b) and the fact that all Jordan blocks of a diagonalizable matrix are of size 1.

Corollary 3.1

Let ${\varvec{I}} - {\varvec{A}}$ is a diagonalizable matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag(\epsilon _{1}, \epsilon _{2}, \dots , \epsilon _{d}) {\varvec{P^{-1}}}$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. If $\epsilon _{i}$’s are real and positive, then

  $$\begin{aligned} \bar{F}_{L}(t) = \left( \frac{\alpha ^{k-\nu }}{\Gamma _{\nu }(k,\alpha l)} \right) {\varvec{a}} {\varvec{P}} {\varvec{D(t)}} {\varvec{P^{-1}}} {\varvec{e}}, \end{aligned}$$

  where

  $$\begin{aligned}{} & {} {\varvec{D(t)}} = diag\left( \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{1} + \alpha )l)}{(\lambda t \epsilon _{1} + \alpha )^{k-\nu }}, \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{2} + \alpha )l)}{(\lambda t \epsilon _{2} + \alpha )^{k-\nu }},\dots , \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{d} + \alpha )l)}{(\lambda t \epsilon _{d} + \alpha )^{k-\nu }} \right) , \\{} & {} \quad \quad t \ge 0. \end{aligned}$$

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,d$, then

  $$\begin{aligned} \bar{F}_{L}(t) = {\varvec{a}} {\varvec{P}} {\varvec{{\widetilde{D}}(t)}} {\varvec{P^{-1}}} {\varvec{e}}, \end{aligned}$$

  where

  $$\begin{aligned} {\varvec{{\widetilde{D}}(t)}} = diag\left( \left( \frac{b}{t \epsilon _{1}+ b }\right) ^{\beta }, \left( \frac{b}{t \epsilon _{2} + b }\right) ^{\beta }, \dots , \left( \frac{b}{t \epsilon _{d} + b }\right) ^{\beta } \right) , \quad t \ge 0. \end{aligned}$$

  Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,d$, then the above result holds for all $\beta \in {\mathbb {R}}$.$\Box$

In the next theorem, we give an alternate representation of the survival function of the above corollary. This will be used in proving the results given in the forthcoming subsections. The proof of part (ii) can be done in the same line as in part(i) and hence, omitted.

Theorem 3.2

Let $M \sim DPH({\varvec{a}}, {\varvec{A}})$, where ${\varvec{A}}$ is a matrix of size d. Further, let ${\varvec{I}} - {\varvec{A}}$ be a diagonalizable matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag(\epsilon _{1}, \epsilon _{2}, \dots , \epsilon _{d}) {\varvec{P^{-1}}}$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. If $\epsilon _{i}$’s are real and positive, then

  $$\begin{aligned} \bar{F}_{L}(t) = \sum _{i=1}^{d}c_{i} \left( \frac{\alpha }{\lambda t \epsilon _{i} + \alpha } \right) ^{k - \nu } \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{i} + \alpha )l)}{\Gamma _{\nu }(k,\alpha l) },\quad t \ge 0, \end{aligned}$$

  where $c_{i}$’s are real constants such that $\sum _{i=1}^{d} c_{i} = 1$.

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,d$, then

  $$\begin{aligned} \bar{F}_{L}(t) = \sum _{i=1}^{d} g_{i} \left( \frac{b}{t \epsilon _{i}+ b }\right) ^{\beta }, \quad t \ge 0, \end{aligned}$$

  where $g_{i}$’s are complex constants such that $\sum _{i=1}^{d} g_{i} = 1$. Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,d$, then the above result holds for all $\beta \in {\mathbb {R}}$,; here $g_{i}$’s are real constants.

Remark 3.1

In Theorem 3.2 (i), it is given that $\sum _{i=1}^{d} c_{i} = 1$. Thus, if $c_{i} \ge 0$, for all i, then L is a mixture of d random variables $V_{1}, V_{2},\dots , V_{d}$ with reliability functions given by

$$\begin{aligned} \bar{F}_{V_{i}}(t) = \left( \frac{\alpha }{\lambda t \epsilon _{i} + \alpha } \right) ^{k - \nu } \frac{\Gamma _{\nu }( k, (\lambda t \epsilon _{i} + \alpha )l)}{\Gamma _{\nu }(k,\alpha l) }, \quad i = 1,2,\dots ,d. \end{aligned}$$

Similarly, in Theorem 3.2 (ii), if $g_{i} \ge 0$, for all i, then L is a mixture of d Pareto distributed random variables with the sets of parameters $\{\beta , b/\epsilon _{i}\}$, $i = 1,2,\dots ,d$. Note that, the Pareto distribution has the DFR (decreasing failure rate) property. Thus, if $g_{i} \ge 0$, for all i, then L is also DFR. In particular, if ${\varvec{I}} - {\varvec{A}}$ is a diagonal matrix, then $g_{i} \ge 0$, for all i, and hence, L is DFR in this case.

3.2 Failure rate function

The proof of part (ii) of the following theorem, where the failure rate $r_{L}(t)$ for a system subject to the defined shock process is derived, can be performed in the same line as in part (i) and hence, omitted.

Theorem 3.3

Let $M \sim DPH({\varvec{a}}, {\varvec{A}})$, where ${\varvec{A}}$ is a square matrix of size d. Further, let $\epsilon _{1},\epsilon _{2},\dots ,\epsilon _{r}$ be eigenvalues of the matrix ${\varvec{I}} - {\varvec{A}}$ with the algebraic multiplicity $d_{1},d_2,\dots ,d_r$, respectively, and let ${\varvec{P}}$ be a non-singular matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag({\varvec{J_{1}}}, {\varvec{J_{2}}}, \dots , {\varvec{J_{r}}}) {\varvec{P^{-1}}}$, where ${\varvec{J_{i}}}$ is the $\epsilon _{i}$-Jordan block with dimension $d_{i}$, $i = 1,2,\dots ,r$, and $\sum _{i=1}^{r}d_{i} = d$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. If every eigenvalue of the matrix ${\varvec{I}} - {\varvec{A}}$ is real and positive, then

  $$\begin{aligned} r_{L}(t) = \lambda \left( \frac{{\varvec{a}} {\varvec{P}} diag({\varvec{H_{1}(t)}}, {\varvec{H_{2}(t)}}, \dots , {\varvec{H_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} diag({\varvec{D_{1}(t)}}, {\varvec{D_{2}(t)}}, \dots , {\varvec{D_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}} \right) , \end{aligned}$$

  where ${\varvec{D_{i}(t)}}$ is same as in Theorem 3.1, and ${\varvec{H_{i}(t)}} = {\varvec{E_{i}(t)}} - {\varvec{S_{i}(t)}}$,

  $$\begin{aligned}{} & {} {\varvec{E_{i}(t)}} = \begin{pmatrix} \frac{\epsilon _{i} \Gamma _{\nu }( k +1, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k + 1-\nu }} &{} \frac{ - \lambda t \epsilon _{i} \Gamma _{\nu }(k + 2, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k + 2-\nu }} &{} \dots &{} \frac{(-\lambda t)^{d_{i}-1} \epsilon _{i} \Gamma _{\nu }(k + d_{i}, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-1)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i}-\nu }} \\ 0 &{} \frac{\epsilon _{i} \Gamma _{\nu }( k + 1, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k + 1-\nu }} &{} \dots &{} \frac{(-\lambda t)^{d_{i}-2} \epsilon _{i} \Gamma _{\nu }(k + d_{i} -1, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-2)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i} -1-\nu }}\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \frac{\epsilon _{i} \Gamma _{\nu }( k + 1, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k+ 1-\nu }} \end{pmatrix}_{d_{i}\times d_{i}}, \\{} & {} \quad {\varvec{S_{i}(t)}} = \begin{pmatrix} 0 &{} \frac{ - \Gamma _{\nu }(k + 1, (\lambda t \epsilon _{i} + \alpha )l)}{(\lambda t \epsilon _{i} + \alpha )^{k + 1-\nu }} &{} \dots &{} \frac{-(-\lambda t)^{d_{i}-2} \Gamma _{\nu }(k + d_{i} -1, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-2)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i} -1-\nu }} \\ 0 &{} 0 &{} \dots &{} \frac{-(-\lambda t)^{d_{i}-3} \Gamma _{\nu }(k + d_{i} -2, (\lambda t \epsilon _{i} + \alpha )l)}{(d_{i}-3)! (\lambda t \epsilon _{i} + \alpha )^{k+d_{i} -2-\nu }}\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} 0 \end{pmatrix}_{d_{i}\times d_{i}}, \end{aligned}$$

  $t \ge 0$, $i = 1,2,\dots ,r$.

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,r$, then

  $$\begin{aligned} r_{L}(t) = \frac{1}{b} \left( \frac{{\varvec{a}} {\varvec{P}} diag({\varvec{{\widetilde{H}}_{1}(t)}}, {\varvec{{\widetilde{H}}_{2}(t)}}, \dots , {\varvec{{\widetilde{H}}_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}}{\varvec{P}} diag({\varvec{{\widetilde{D}}_{1}(t)}}, {\varvec{{\widetilde{D}}_{2}(t)}}, \dots , {\varvec{{\widetilde{D}}_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}} \right) , \end{aligned}$$

  where ${\varvec{{\widetilde{D}}_{i}(t)}}$ is the same as in Theorem 3.1, and ${\varvec{{\widetilde{H}}_{i}(t)}} = {\varvec{{\widetilde{E}}_{i}(t)}} - {\varvec{{\widetilde{S}}_{i}(t)}}$,

  $$\begin{aligned}{} & {} {\varvec{\widetilde{E_{i}}(t)}} = \begin{pmatrix} \epsilon _{i} \Gamma (\beta + 1)\left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1} &{} - \epsilon _{i} \Gamma (\beta + 2) \left( \frac{t}{t \epsilon _{i} + b}\right) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1}&{} \dots &{} \frac{ \epsilon _{i} \Gamma (\beta + d_{i})}{(d_{i}-1)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-1} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1} \\ 0 &{} \epsilon _{i} \Gamma (\beta + 1) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1} &{} \dots &{} \frac{ \epsilon _{i} \Gamma (\beta + d_{i}-1)}{(d_{i}-2)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-2} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1}\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \epsilon _{i} \Gamma (\beta + 1) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1} \end{pmatrix}_{d_{i} \times d_{i}}, \\{} & {} \quad {\varvec{\widetilde{S_{i}}(t)}} = \begin{pmatrix} 0 &{} - \Gamma (\beta + 1) \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1}&{} \dots &{} \frac{- \Gamma (\beta + d_{i}-1)}{(d_{i}-2)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-2} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1} \\ 0 &{} 0 &{} \dots &{} \frac{- \Gamma (\beta + d_{i}-2)}{(d_{i}-3)!} \left( \frac{- t}{t \epsilon _{i} + b}\right) ^{d_{i}-3} \left( \frac{b}{t \epsilon _{i} + b}\right) ^{\beta + 1}\\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} 0 \end{pmatrix}_{d_{i} \times d_{i}}, \end{aligned}$$

  $t \ge 0$,$i = 1,2,\dots ,r$. Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,r$, then the above result holds for all $\beta \in {\mathbb {R}}$.$\Box$

The following corollary immediately follows from Theorem 3.3, Remark 2.2(b) and the fact that all Jordan blocks of a diagonalizable matrix are of size 1.

Corollary 3.2

Let ${\varvec{I}} - {\varvec{A}}$ is a diagonalizable matrix such that ${\varvec{I}} -{\varvec{A}} = {\varvec{P}} diag (\epsilon _{1}, \epsilon _{2}, \dots , \epsilon _{d}) {\varvec{P^{-1}}}$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. If $\epsilon _{i}$’s are real and positive, then

  $$\begin{aligned} r_{L}(t) = \lambda \left( \frac{ {\varvec{a}} {\varvec{P}} {\varvec{H(t)}} {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} {\varvec{D(t)}} {\varvec{P^{-1}}} {\varvec{e}}}\right) , \end{aligned}$$

  where ${\varvec{D(t)}}$ is the same as in Corollary 3.1 and

  $$\begin{aligned} {\varvec{H(t)}}= & {} diag\left( \frac{ \epsilon _{1} \Gamma _{\nu }( k +1, (\lambda t \epsilon _{1} + \alpha )l)}{(\lambda t \epsilon _{1} + \alpha )^{k + 1-\nu }}, \frac{\epsilon _{2} \Gamma _{\nu }( k +1, (\lambda t \epsilon _{2} + \alpha )l)}{(\lambda t \epsilon _{2} + \alpha )^{k + 1-\nu }},\dots ,\right. \\{} & {} \qquad \left. \frac{\epsilon _{d} \Gamma _{\nu }( k +1, (\lambda t \epsilon _{d} + \alpha )l)}{(\lambda t \epsilon _{d} + \alpha )^{k + 1-\nu }} \right) , \end{aligned}$$

  for $t \ge 0$.

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,d$, then

  $$\begin{aligned} r_{L}(t) = \frac{\beta }{b} \left( \frac{ {\varvec{a}} {\varvec{P}} {\varvec{{\widetilde{H}}(t)}} {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} {\varvec{{\widetilde{D}}(t)}} {\varvec{P^{-1}}} {\varvec{e}}}\right) , \end{aligned}$$

  where ${\varvec{{\widetilde{D}}(t)}}$ is the same as in Corollary 3.1 and

  $$\begin{aligned}{} & {} {\varvec{{\widetilde{H}}(t)}} = diag\left( \epsilon _{1} \left( \frac{b}{t \epsilon _{1} + b }\right) ^{\beta + 1}, \epsilon _{2} \left( \frac{b}{t \epsilon _{2}+ b}\right) ^{\beta + 1}, \dots , \epsilon _{d} \left( \frac{b}{t \epsilon _{d} + b}\right) ^{\beta + 1} \right) ,\\{} & {} \quad \quad t \ge 0. \end{aligned}$$

  Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,d$, then the above result holds for all $\beta \in {\mathbb {R}}$.

     From Remark  3.1, we already know that the cdf $F_L(t)$ has a decreasing failure rate. The following intuitively clear theorem (as the failure rate of the Pareto distribution asymptotically tends to 0) shows that the failure rate that corresponds to $F_L(t)$ also tends to 0.

Theorem 3.4

Let ${\varvec{M}} \sim DPH({\varvec{a}}, {\varvec{A}})$, where ${\varvec{A}}$ is a square matrix of size d. Further, let ${\varvec{I}} - {\varvec{A}}$ be a diagonalizable matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag(\epsilon _{1}, \epsilon _{2}, \dots , \epsilon _{d}) {\varvec{P^{-1}}}$. Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,d$, then

$$\begin{aligned} \lim _{t \rightarrow \infty } r_{L}(t) = 0. \end{aligned}$$

Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,d$, then this result holds for all $\beta \in {\mathbb {R}}$.

3.3 Mean residual lifetime and mean lifetime

The mean residual life function is also an important reliability characteristic for a system operating under shocks. In the following theorem, we derive the corresponding expressions for the defined shock model. We prove only part (i), whereas part (ii) can be performed in a similar way and, hence, omitted.

Theorem 3.5

Let $M \sim DPH({\varvec{a}}, {\varvec{A}})$, where ${\varvec{A}}$ is a square matrix of size d. Further, let $\epsilon _{1},\epsilon _{2},\dots ,\epsilon _{r}$ be eigenvalues of the matrix ${\varvec{I}} - {\varvec{A}}$ with the algebraic multiplicity $d_{1},d_2,\dots ,d_r$, respectively, and let ${\varvec{P}}$ be a non-singular matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag({\varvec{J_{1}}}, {\varvec{J_{2}}}, \dots , {\varvec{J_{r}}}) {\varvec{P^{-1}}}$, where ${\varvec{J_{i}}}$ is the $\epsilon _{i}$-Jordan block with dimension $d_{i}$, $i = 1,2,\dots ,r$, and $\sum _{i=1}^{r}d_{i} = d$.

• (i) Assume that shocks occur according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$, $k > 1$. If every eigenvalue of the matrix ${\varvec{I}} - {\varvec{A}}$ is real and positive, then

  $$\begin{aligned} E(L-t|L>t) = \frac{ {\varvec{a}} {\varvec{P}} diag({\varvec{R_{1}(t)}}, {\varvec{R_{2}(t)}}, \dots , {\varvec{R_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} diag({\varvec{D_{1}(t)}}, {\varvec{D_{2}(t)}}, \dots , {\varvec{D_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}, \end{aligned}$$

  where ${\varvec{D_{i}(t)}}$ is the same as in Theorem 3.1, and

  $$\begin{aligned} {\varvec{R_{i}(t)}}= & {} \begin{pmatrix} \phi (0,t,\epsilon _{i}) &{}-\phi (1,t,\epsilon _{i}) &{} \dots &{} \frac{(-1)^{d_{i}-1} \phi (d_{i}-1,t,\epsilon _{i})}{(d_{i}-1)!} \\ 0 &{} \phi (0,t,\epsilon _{i}) &{} \dots &{} \frac{(-1)^{d_{i}-2} \phi (d_{i}-2,t,\epsilon _{i})}{(d_{i}-2)!} \\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \phi (0,t,\epsilon _{i}) \end{pmatrix}_{d_{i}\times d_{i}}, \\ \phi (n,t,\epsilon _{i})= & {} \frac{ (\lambda t)^{n} \Gamma _{\nu }(k+n-1, (\alpha + \lambda t \epsilon _{i})l)}{\lambda \epsilon _{i} (\alpha + \lambda t \epsilon _{i})^{n+k-\nu -1}} \\{} & {} + \frac{(\lambda t)^{n}}{\lambda \epsilon _{i}} \sum _{i=1}^{n} \frac{n (n-1) \dots (n-i+1)}{(\lambda t \epsilon _{i})^{i}} \frac{\Gamma _{\nu }(n+k-i-1,(\alpha +\lambda t \epsilon _{i})l)}{(\alpha +\lambda t \epsilon _{i})^{n+k-\nu -i-1}}, \end{aligned}$$

  for $t \ge 0$, $i = 1,2,\dots ,r$ and $n \in {\mathbb {N}} \cup \{0\}$.

• (ii) Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$, $\beta > 1$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,r$, then

  $$\begin{aligned} E(L-t|L>t) = \frac{ {\varvec{a}} {\varvec{P}} diag({\varvec{{\widetilde{R}}_{1}(t)}}, {\varvec{{\widetilde{R}}_{2}(t)}}, \dots , {\varvec{{\widetilde{R}}_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} diag({\varvec{{\widetilde{D}}_{1}(t)}}, {\varvec{{\widetilde{D}}_{2}(t)}}, \dots , {\varvec{{\widetilde{D}}_{r}(t)}}) {\varvec{P^{-1}}} {\varvec{e}}}, \end{aligned}$$

  where ${\varvec{{\widetilde{D}}_{i}(t)}}$ is the same as in Theorem 3.1, and

  $$\begin{aligned} {\varvec{\widetilde{R_{i}}(t)}}= & {} \begin{pmatrix} \psi (0,t,\epsilon _{i}) &{}-\psi (1,t,\epsilon _{i}) &{} \dots &{} \frac{(-1)^{d_{i}-1} \psi (d_{i}-1,t,\epsilon _{i})}{(d_{i}-1)!} \\ 0 &{} \psi (0,t,\epsilon _{i}) &{} \dots &{} \frac{(-1)^{d_{i}-2} \psi (d_{i}-2,t,\epsilon _{i})}{(d_{i}-2)!} \\ \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \psi (0,t,\epsilon _{i}) \end{pmatrix}, \\{} & {} \quad \psi (n,t,\epsilon _{i}) = (b/\epsilon _{i}) \Gamma (\beta + n -1) \left( \frac{t}{b+ t\epsilon _{i}}\right) ^{n} \left( \frac{b}{b+t\epsilon _{i}}\right) ^{\beta -1} \\{} & {} \quad + (b/\epsilon _{i}) \left( \frac{t}{b + t\epsilon _{i}}\right) ^{n} \sum _{i=1}^{n} \frac{b^{i} n (n-1) \dots (n-i+1)}{( t \epsilon _{i})^{i} \Gamma (\beta + n - i - 1)} \left( \frac{b}{b + t \epsilon _{i}}\right) ^{\beta -i-1}, \end{aligned}$$

  for $t \ge 0$, $i = 1,2,\dots ,r$, and $n \in {\mathbb {N}} \cup \{0\}$. Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,r$, then the above result holds for all $\beta \in {\mathbb {R}}$.$\Box$

The following corollary immediately follows from Theorem 3.5, Remark 2.2(b) and the fact that all Jordan blocks of a diagonalizable matrix are of size 1.

Corollary 3.3

Let ${\varvec{I}} - {\varvec{A}}$ be a diagonalizable matrix such that ${\varvec{I}} - {\varvec{A}} = {\varvec{P}} diag(\epsilon _{1}, \epsilon _{2}, \dots , \epsilon _{d}) {\varvec{P^{-1}}}$.

• (i) Assume shocks occur on the system according to the HPGGP with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$, $k > 1$. If $\epsilon _{i}$’s are real and positive, then

  $$\begin{aligned} E(L-t|L>t) = \frac{ {\varvec{a}} {\varvec{P}} {\varvec{R(t)}} {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} {\varvec{D(t)}} {\varvec{P^{-1}}} {\varvec{e}}}, \end{aligned}$$

  where ${\varvec{D(t)}}$ is the same as in Corollary 3.1, and

  $$\begin{aligned}{} & {} {\varvec{R(t)}} = diag( \phi (0,t,\epsilon _{i}), \phi (0,t,\epsilon _{2}), \dots , \phi (0,t,\epsilon _{d})), \\{} & {} \phi (0,t,\epsilon _{i}) = \frac{ \Gamma _{\nu }(k-1, (\alpha + \lambda t \epsilon _{i})l)}{\lambda \epsilon _{i} (\alpha + \lambda t \epsilon _{i})^{k-\nu -1}}, \quad t \ge 0, i = 1,2,\dots ,d. \end{aligned}$$

• (ii) Assume that shocks occur on the system according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta \in {\mathbb {N}}$, $\beta > 1$. If $Re(\epsilon _{i}) > 0$, for all $i = 1,2,\dots ,d$, then

  $$\begin{aligned} E(L-t|L>t) = \frac{ {\varvec{a}} {\varvec{P}} {\varvec{{\widetilde{R}}(t)}} {\varvec{P^{-1}}} {\varvec{e}}}{ {\varvec{a}} {\varvec{P}} {\varvec{{\widetilde{D}}(t)}} {\varvec{P^{-1}}} {\varvec{e}}}, \end{aligned}$$

  where ${\varvec{{\widetilde{D}}(t)}}$ is the same as in Corollary 3.1, and

  $$\begin{aligned} {\varvec{{\widetilde{R}}(t)}}= & {} diag\left( (b/\epsilon _{i}) \left( \frac{1}{\beta -1}\right) \left( \frac{b}{b+t\epsilon _{i}}\right) ^{\beta -1},\dots ,\right. \\{} & {} \quad \left. (b/\epsilon _{i}) \left( \frac{1}{\beta -1}\right) \left( \frac{b}{b+t\epsilon _{i}}\right) ^{\beta -1}\right) , \end{aligned}$$

  for $t \ge 0$. Moreover, if $Im(\epsilon _{i}) = 0$, for all $i = 1,2,\dots ,d$, then the above result holds for all $\beta \in {\mathbb {R}}$.$\Box$

In the following theorem, the mean lifetime of the system is given for the defined shock model. Although it is a specific case for the mean residual lifetime, it is more efficient and simple to obtain this result directly and under weaker assumptions.

Theorem 3.6

Let $M \sim DPH({\varvec{a}},{\varvec{A}})$. Assume that shocks occur according to the HPGGP with parameters $\{\lambda , \nu ,k,\alpha ,l\}$, $k > 1$. Then,

$$\begin{aligned} E(L) = \frac{\alpha ({\varvec{a}} ({\varvec{I}} - {\varvec{A}})^{-1} {\varvec{e}})}{\lambda } \frac{ \Gamma _{\nu }(k-1, \alpha l)}{\Gamma _{\nu }(k, \alpha l)}. \end{aligned}$$

Corollary 3.4

Assume that shocks occur according to the Pólya process with the set of parameters $\{\beta , b\}$, $\beta > 1$. Then,

$$\begin{aligned} E(L) = \frac{b ({\varvec{a}} ({\varvec{I}} - {\varvec{A}})^{-1} {\varvec{e}})}{\beta - 1}. \end{aligned}$$

3.4 Stochastic comparisons

In this subsection, we study some stochastic comparisons for systems for the defined shock model. First, we recall the following definition

Definition 3.2

A random variable V is said to be smaller than another random variable W in the usual stochastic order, denoted by $V \le _{st} W$, if $\bar{F}_{V}(t) \le \bar{F}_{W}(t)$ for all $t >0$. $\Box$

Theorem 3.7

Let $L_{1} = \Sigma _{i=1}^{M_{1}} X_{i}$ and $L_{2} = \Sigma _{i=1}^{M_{2}} X_{i}$ be the lifetimes of two systems subject to the same HPGGP of shocks with the set of parameters $\{\lambda , \nu ,k,\alpha ,l\}$. Further, let $M_{i} \sim DPH({\varvec{\gamma _{i}}}, {\varvec{A_{i}}})$, where ${\varvec{I}} - {\varvec{A_{i}}}$ is a diagonal matrix of order d with positive and real eigenvalues $\epsilon ^{(i)}_{1}, \epsilon ^{(i)}_{2}, \dots , \epsilon ^{(i)}_{d}$ such that $\epsilon ^{(i)}_{1} \le \epsilon ^{(i)}_{2} \le \dots \le \epsilon ^{(i)}_{d}$, $i=1,2$. If ${\varvec{\gamma _{1}}} \ge _{st} {\varvec{\gamma _{2}}}$ and $\epsilon ^{(1)}_{j} \ge \epsilon ^{(2)}_{j}$, for all $j = 1,2,\dots ,d$.

Then, $L_{1} \le _{st} L_{2}.$ $\Box$

The following corollary immediately follows from Theorem 3.7.

Corollary 3.5

The following results hold true.

• (a) Assume that ${\varvec{A_{1}}} = {\varvec{A_{2}}}$. If ${\varvec{\gamma _{1}}} \ge _{st} {\varvec{\gamma _{2}}}$, then $L_{1} \le _{st} L_{2}$;

• (b) Assume that ${\varvec{\gamma _{1}}} = {\varvec{\gamma _{2}}}$. If $\epsilon ^{(1)}_{j} \ge \epsilon ^{(2)}_{j}$, for all $j = 1,2,\dots ,d$, then $L_{1} \le _{st} L_{2}$;

• (c) Assume that shocks occur on the system according to the Pólya process with the set of parameters $\{\beta , b\}$. If ${\varvec{\gamma _{1}}} \ge _{st} {\varvec{\gamma _{2}}}$ and $\epsilon ^{(1)}_{j} \ge \epsilon ^{(2)}_{j}$, for all $j = 1,2,\dots ,d$, then $L_{1} \le _{st} L_{2}.$

Theorem 3.8

Let $L_{1} = \Sigma _{i=1}^{M} X^{(1)}_{i}$ and $L_{2} = \Sigma _{i=1}^{M} X^{(2)}_{i}$ be the lifetimes of two systems subject to different HPGGP shock processes with parameter sets $\{\lambda _{1}, \nu _{1},k_{1},\alpha _{1},l_{1}\}$ and $\{\lambda _{2}, \nu _{2},k_{2},\alpha _{2},l_{2}\}$, respectively, where $X^{(1)}_{i}$ and $X^{(2)}_{i}$ are the ith inter-arrival times of shock processes that impact the first and the second systems, respectively. Further, let $M \sim DPH({\varvec{a}}, {\varvec{A}})$. Suppose that $\lambda _{1} \le \lambda _{2}$, and one of the following conditions holds.

• (a) $\alpha _{1} = \alpha _{2}$, $k_{2} \ge k_{1}$ and $k_{2} - k_{1} \ge \nu _{2} - \nu _{1}$;

• (b) $\alpha _{2} < \alpha _{1}$ and $(\alpha _{1} - \alpha _{2} + k_{2} - k_{1} + \nu _{1} - \nu _{2})^{2} - 4 (\alpha _{1} - \alpha _{2})(k_{2} - k_{1}) \le 0$;

• (c) $\alpha _{2} < \alpha _{1}$ and $\alpha _{1} - \alpha _{2} + k_{2} - k_{1} + \nu _{1} - \nu _{2} \ge 0$.

Then, $L_{2} \le _{st} L_{1}$. $\Box$

The following corollary immediately follows from Theorem 3.8 and Remark 2.2(b).

Corollary 3.6

Assume that shocks occur on the first and the second systems according to the Pólya processes with parameter sets $\{\beta _{1}, b_{1}\}$ and $\{\beta _{2}, b_{2}\}$, respectively. If $b_{1} \ge b_{2}$ and $\beta _{1} \le \beta _{2}$, then $L_{2} \le _{st} L_{1}$.

4 Illustrative examples

In this section, we discuss some numerical examples to illustrate the results given in Sect. 3.

Example 4.1

Consider the generalized extreme shock model I as mentioned in Subsect. 2.3. We assume that $n = 3$, $\theta _{1} = 0.6$, $\theta _{2} = 0.3$ and $\theta _{3} = 0.1$. In Fig. 1a, c and e, we plot the survival function, the failure rate function and the mean residual life function, respectively, for different $p_i$’s and for fixed process parameters $b = 2$ and $\beta = 4$. Further, in Fig. 1b, d, f we plot the survival function, the failure rate function and the mean residual life function, respectively, for different process parameters and for fixed $p_{1} = 0.1$, $p_{2} = 0.2$ and $p_{3} = 0.3$. Clearly, Fig. 1a and b demonstrates the results given in Corollaries 3.5 and 3.6, respectively. Figure 1c and d demonstrates the result given in Theorem 3.4. Moreover, from these figures, we see that the failure rate of the system is DFR (see Remark 3.1). Further, Fig. 1e and f shows that the system’s mean residual life function is nondecreasing.

We plot Fig. 2 that illustrates the impact of the dependent increments property of the HPGGP on different reliability measures derived in the previous sections. In Fig. 2a, we compare systems’ survival functions for the HPP (the process with independent increments) and the HPGGP (the process with dependent increments) shock processes. We clearly see that survival functions for dependent and independent increments cases are significantly different. Similarly, in Fig. 2b and c, we compare systems’ failure rate functions and mean residual lifetime functions, respectively, for the HPP and the HPGGP cases. These figures also show a similar conclusion.

Fig. 1

Plot of the system’s survival, failure rate and mean residual life functions for the generalized extreme shock model I

Fig. 2

Plot of the system’s survival, failure rate and mean residual life functions for HPP and HPGGP shock processes

Example 4.2

Consider the classical run shock model (with $r = 2$) as mentioned in Subsect. 2.3. Then,

$$\begin{aligned} {\varvec{a_{rm}}} = (1,0), \quad {\varvec{A_{rm}}} = \begin{pmatrix} 1 - p &{}\quad p \\ 1 - p &{}\quad 0 \end{pmatrix}. \end{aligned}$$

Clearly, ${\varvec{I}} - {\varvec{A_{rm}}}$ is diagonalizable and has positive and real eigenvalues for all $p \in [0,1)$. For $p = 0.2$, eigenvalues of the matrix ${\varvec{I}} -{\varvec{ A_{rm}}}$ are given by $\epsilon ^{(1)}_{1} = 1.16568542$ and $\epsilon ^{(1)}_{2} = 0.03431458$. Similarly, for $p = 0.3$, eigenvalues of the matrix ${\varvec{I}} - {\varvec{A_{rm}}}$ are given by $\epsilon ^{(2)}_{1} = 1.22662813$ and $\epsilon ^{(2)}_{2} = 0.07337187$. Clearly, all eigenvalues, for $p = 0.2$ and $p = 0.3$, are real and positive. In Fig. 3a,  c and  e, we plot the system’s survival, failure rate and mean residual life functions, respectively, for different p values by assuming that shock process is the HPGGP with the set of parameters $\{2,2,2,2,3\}$. In Fig. 3b, we plot the system’s survival function for fixed $p = 0.2$ and for the HPGGP with different values of parameters $\{\lambda , \nu , k, \alpha , l\}$. In Fig. 3d and  f, we plot the system’s failure rate and mean residual life functions for fixed $p = 0.3$ and for the HPGGP with different values of parameters $\{\lambda , \nu , k, \alpha , l\}$. Since $\epsilon ^{(1)}_{j} \le \epsilon ^{(2)}_{j}$ for $j = 1,2$, Fig. 3a illustrates the result given in Theorem 3.7. Further, Fig. 3b illustrates the result given in Theorem 3.8 as condition (a) of Theorem 3.8 holds. Figure 3c and  d demonstrates that the failure rate of the system asymptotically tends to 0. Figure 3e and  f shows that the system’s mean residual life function is nondecreasing.

Fig. 3

Plot of the system’s survival, failure rate, and mean residual life functions for the run shock model ($r = 2$)

Example 4.3

Consider the classical run shock model (with $r = 3$) as mentioned in Subsection 2.3. Then,

$$\begin{aligned} {\varvec{a_{rm}}} = (1,0,0), \quad {\varvec{A_{rm}}} = \begin{pmatrix} 1-p &{}\quad p &{}\quad 0 \\ 1-p &{} \quad 0 &{}\quad p \\ 1-p &{} \quad 0 &{}\quad 0 \end{pmatrix}. \end{aligned}$$

We assume $p=0.2$ and 0.3. For $p=0.2$, eigenvalues of the matrix ${\varvec{I}} - {\varvec{A_{rm}}}$ are $1.096736514+0.15117 i$, $1.096736514-0.15117i$ and 0.006526972, where $i=\sqrt{-1}$. Similarly, for $p = 0.3$, eigenvalues of the matrix ${\varvec{I}} - {\varvec{A_{rm}}}$ are $1.13995689+0.211432 i$, $1.13995689-0.211432 i$ and 0.02008622. Clearly, all eigenvalues, for both $p = 0.2$ and $p = 0.3$, are either real or complex with positive real parts. In Fig. 4a, c and e, we plot the system’s survival, failure rate and mean residual life functions, respectively, for different p values and for the Pólya shock process with parameters $b = 2$ and $\beta = 4$. In Fig. 4b, we plot the system’s survival function for fixed $p = 0.2$ and for the Pólya shock process with different values of parameters b and $\beta$. In Fig. 4d and f, we plot the system’s failure rate and mean residual life functions for fixed $p = 0.4$ and for the Pólya shock process with different values of parameters b and $\beta$. Figure 4a shows that an increment in the parameter p decreases the system’s survivability. Further, Fig. 4b illustrates the result given in Corollary 3.6. Figure 4c and d demonstrates the result given in Theorem 3.4. Further, Fig. 4e shows that an increment in the parameter p decreases the system’s mean residual life.

Fig. 4

Plot of the system’s survival, failure rate, and mean residual life functions for the run shock model $(r = 3)$

Example 4.4

Consider the generalized run shock model I with $k_{1} = 2$ and $k_{2} = 1$. Then,

$$\begin{aligned} {\varvec{a_{grm1}}} = (1,0)\text { and } {\varvec{A_{grm1}}} = \begin{pmatrix} p_{1} &{}\quad p_{2} \\ p_{1} &{} \quad 0 \end{pmatrix}. \end{aligned}$$

We assume $p_{1} = 0.7$ and $p_{2} = 0.1$. In Fig. 5a, we plot the system’s survival function for the Pólya shock process with different parameter values. This illustrates the result given in Corollary 3.6. In Fig. 5b and c, we plot the system’s failure rate and mean residual life functions for fixed $p_{1} = 0.6$ and $p_{2} = 0.3$. Figure 5b demonstrates the result given in Theorem 3.4. From Fig. 5c, we can see that mean residual life function is nondecreasing.

Fig. 5

Plot of system’s survival, failure rate, and mean residual life functions for generalized run shock model I

5 Application: optimal mission duration

In this section, we consider the optimal mission duration policy (see Finkelstein and Levitin (2018)) for a system subject to shocks, whereas its failure is defined in accordance with the classical run shock model. The same can be done for other models as well. To avoid a failure of a system during a mission of a fixed duration t (that can cause substantial losses), in many real-life scenarios, it may be a reasonable strategy to abort the mission before its completion. For example, consider a chemical production system that has to supply a contracted amount of some commodity. If the contract is fulfilled (i.e., the predetermined amount of the commodity is supplied), the producer receives an award. However, external impacts can cause a failure of the system during the operation. The failure can destroy the system (in some cases, it also destroys the commodity produced so far). If the producer estimates the risk associated with the failure and decides to terminate the operation, it can pay certain penalties for default of the contract. Depending on the nature of the commodity, the producer can or cannot sell the amount of commodity produced before the mission termination (see Finkelstein and Levitin (2018)). The mission abort usually results in a reward that depends on the system’s operation time and the corresponding penalty. On the other hand, the mission completion results in an additional reward. Moreover, the failure of the system during the mission also results in a penalty because it incurs additional costs due to failure of a mission. The decision about the mission termination at any time $\tau$ should be made if the profit in case of termination exceeds the expected profit in case of its continuation with respect to risk associated with the system future failure (see Finkelstein and Levitin (2018)).

Assumptions:

• 1. Let a new system be installed at time $t=0$, and let its lifetime be modeled by the classical run shock model. Further, let L and $t \; ( >0)$ denote the lifetime of the system and the mission duration, respectively. The mission can be terminated at any time $\tau \in (0,t]$.

• 2. Shocks occur according to the Pólya process with parameter set $\{\beta , b\}$.

• 3. The profit C(t) is obtained when the mission is completed (i.e., the system does not fail during the mission or the mission is not aborted in [0, t])). The per time unit reward, when the system is operating, is $c_{p}$ and the per time unit operational cost is $c_{0}$, where $c_{0} < c_{p}$.

• 4. The penalty $C_{f}$ is imposed if the system fails during the mission. In case of the premature termination, the fixed penalty $C_{pt}$ ($C_{pt} < C_{f}$) is administrated. Further, $C_{r}$ is an additional reward for the mission completion.

• 5. The reward after the failure is discarded.

Based on the aforementioned assumptions, the profit C(t) for the mission completion is given by (see Finkelstein and Levitin (2018))

$$\begin{aligned} C(t) = (c_{p} - c_{0})t + C_{r}. \end{aligned}$$

Note that the mission is aborted at time $\tau$ if the total profit at termination exceeds the expected profit in case of mission continuation. The profit at termination at time $\tau$ is equal to $(c_{p} - c_{0}) \tau - C_{pt}$. On the other hand, the expected profit in the case of mission continuation is

$$\begin{aligned} \frac{{\bar{F}}_{L}(t)}{{\bar{F}}_{L}(\tau )} \left( \left( c_{p} - c_{0} \right) t + C_{r} \right) - \left( 1 - \frac{{\bar{F}}_{L}(t)}{{\bar{F}}_{L}(\tau )} \right) C_{f}, \end{aligned}$$

where ${{\bar{F}}_{L}(t)}/{{\bar{F}}_{L}(\tau )}$ is the probability that a system will not fail in the remaining mission time given that it is operable at time $\tau$; here ${\bar{F}}_{L}(t)$ is the same as in Theorem 3.1. Thus, if for some $\tau$, the expression

$$\begin{aligned} A(\tau ) {\mathop {=}\limits ^\mathrm{def.}} \frac{{\bar{F}}_{L}(t)}{{\bar{F}}_{L}(\tau )} \left( \left( c_{p} - c_{0} \right) t + C_{r} \right) - \left( 1 - \frac{{\bar{F}}_{L} (t)}{{\bar{F}}_{L} (\tau )} \right) C_{f} - \left( (c_{p} - c_{0}) \tau - C_{pt} \right) \end{aligned}$$

is nonnegative, then the mission should not be terminated at time $\tau$. Clearly, $A(0) \ge 0$, as there is no need to terminate the mission that had just started. Since the expression of $A(\tau )$ is complicated, it is not analytically possible to find out the values of $\tau$ for which $A(\tau )\ge 0$. Thus, we consider the following numerical example.

Fig. 6

Profit comparison function $A(\tau )$ against $\tau \in [0,10]$

Let us assume $t = 10, c_{p} = 2.5, c_{0} = 0.5, C_{r} = 3, C_{f} = 8$, $C_{pt} = 5$. Further, we assume that the run shock model parameters are $r = 2$ and $p = 0.25$, and the process parameters are $b = 1$ and $\beta = 5$. In Fig. 6, we plot the profit comparison function $A(\tau )$ against $\tau \in [0,10]$. This figure shows that $A(\tau )$ is in U-shape. Further, it takes negative values in $\tau \in [0.89, 6.17]$. This implies that the mission should be aborted just at $\tau = 0.89$. In case the mission cannot be aborted at this time (due, for example, to operational or environmental reasons), it may be aborted at any time in the interval (0.89, 3.68], where at $\tau = 3.68$, the function achieves its minimum. Note that there is no other ’optimal’ time to abort except $\tau = 0.89$ in our setting, however, it is not cost-efficient to abort in (3.68, 6.17] and in (6.17, 10] as well, as the plotted function is increasing in these intervals.

6 Concluding remarks

In this paper, we have introduced a general class of shock models when a system lifetime can be written as a random sum of shocks inter-arrival times. The proposed class contains many popular shock models, namely the classical extreme shock model, the generalized extreme shock model, the classical run shock model, the generalized run shock model, specific types of mixed shock models, etc. We study this model under the assumption that shocks occur according to the homogeneous Poisson generalized gamma process (HPGGP). This process is a generalization of many commonly used in applications counting processes, namely the HPP, the Pólya process, etc. This process possesses the dependent increments property that describe real-life phenomena. Moreover, we assume that the distribution of a fatal shock is a discrete PH distribution. As its distribution function and other statistical measures (e.g., mean, variance, etc.) can be written in the matrix forms, the results obtained for the defined class of shock models are also in the matrix forms that can easily be numerically analyzed using different software packages.

For the defined class of shock models, we have derived the major reliability indices (namely the survival function, the failure rate function, the mean residual lifetime function and the mean lifetime) and have discussed some stochastic comparisons with relevant numerical examples. Finally, as an application, the optimal mission duration policy was studied.

Even though the proposed class of shock models contains many popular shock models, there are some other well-known shock models (namely cumulative shock model, $\delta$-shock model) that are not the members of this class. Therefore, the future study of these models under the HPGGP of shocks can constitute a topic for further research.

Appendix: Supplementary data