A Preventive Replacement Policy for a System Subject to Bivariate Generalized Polya Failure Process

Abstract

Numerous studies on preventive maintenance of minimally repaired systems with statistically independent components have been reported in reliability literature. However, in practice, the repair can be worse-than-minimal and the components of a system can be statistically dependent. The existing literature does not cover this important in-practice setting. Therefore, our paper is the first to deal with these issues by modeling dependence in the bivariate set up when a system consists of two dependent parts. We employ the bivariate generalized Polya process to model the corresponding failure and repair process. Relevant stochastic properties of this process have been obtained in order to propose and further discuss the new optimal bivariate preventive maintenance policy with two decision parameters: age and operational history. Moreover, introducing these two parameters in the considered context is also a new feature of the study. Under the proposed policy, the long-run average cost rate is derived and the optimal replacement policies are investigated. Detailed numerical examples illustrate our findings and show the potential efficiency of the obtained results in practice.

1. Introduction

Preventive maintenance (PM) is usually performed on degrading systems in order to decrease the probabilities of failures during operation that can result in substantial losses. As the cost of the PM is smaller than that of a repair upon failure (also taking into account the additional losses due to failures), the corresponding cost-wise optimization problems can be formulated and solved. In this way, an optimal PM time can be obtained that minimizes, e.g., the corresponding cost rate. Thousands of papers devoted to different PM problems and several books entirely dealing with this important in-practice problem have been published in recent decades (see, e.g., the following influential monographs: [1,2,3]).

There have been numerous studies on PM models, where various counting processes are used to model the corresponding failure/repair process. Until now, most of these studies were focused on univariate counting processes of failure/repair such as the nonhomogeneous Poisson process (NHPP) or the renewal process (see [4,5,6,7,8,9] to name a few). For instance, in [4], it was assumed that the failure process follows the NHPP, which means that the corresponding repair type is a minimal repair. In [5], a system subject to two types of failures (minor and catastrophic failures) and repairs (minimal and perfect repairs) was considered for the maintenance optimization. Thus, from the process point of view, it corresponds to the combination of an NHPP and a renewal process.

The well-known minimal repair assumption holds when the corresponding system is composed of a large number of statistically independent components. Hence, its failure rate (FR) ‘is practically unchanged’ after the replacement of the failed component by a new one. However, in real life, the remaining non-failed components are often affected by the failure of a component in a system because it causes additional stress or damage to them. This eventually results in a worse-than-minimal repair of a system as the states of the non-failed components after the minimal repair of the failed component can be ‘worse’ than just before the failure. Some relevant examples of this situation are as follows ([10,11,12]):

(i)

The failure of a still wire cable in a bridge or in an elevator instantaneously increases the stress on the remaining cables and leads to some damages.

(ii)

For a multi-engine airplane, the failure of an engine during flight instantaneously causes increased stress on the non-failed engines.

(iii)

A failure of a pump in a multi-pump hydraulic control system instantly increases the pressure for each non-failed pump.

Recently, as a generalization of the NHPP, a new counting process (called the generalized Polya process (GPP)) has been defined and applied for modeling the univariate failure/repair processes [13]. It is important to note that, under the GPP model, the repair is worse-than-minimal (GPP repair), which makes this process an effective tool for modeling the corresponding optimal PM policies in this case. Specifically, in [12], two periodic PM models were considered and some properties of the optimal policies assuming the GPP repair process were studied. Furthermore, Ref. [14] proposed a generalized replacement policy that already considers the operational history of a system.

The forgoing applies to univariate counting processes, whereas stochastically dependent multivariate series of events arise in many contexts. For some examples in reliability applications, finance, and economics, see [15,16,17,18]. Ref. [17] suggested a general theoretical framework for the multivariate counting process. Recently, new classes of multivariate counting processes have been developed in the literature (see [19,20,21]). Specifically, in [20], the multivariate generalized Polya process (MVGPP) with mathematically tractable properties was defined.

However, to the best of our knowledge, applications of multivariate processes to maintenance models has not yet been developed in the literature. Therefore, in this paper, assuming that the failure process follows the bivariate generalized Polya process (BVGPP) developed in [20], we propose and discuss a new bivariate preventive maintenance policy based on two ‘parameters’: age and operational history. As in [14], where the failure process was univariate, we show the superiority of the proposed policy compared with the original age-based replacement policy for the BVGPP.

In accordance with the foregoing discussion, we want to concisely emphasize the motivation and the novelty of our study:

-

Motivation: Most systems in real life have dependent components, whereas the existing literature does not cover this aspect. Moreover, to the best of our knowledge, until now there have been no studies that consider the PM models with worse-than-minimal repair in multicomponent systems (that also often occurs in practice).

-

Novelty: We employ the bivariate generalized Polya process to model the corresponding failure and repair process. This modeling approach was not considered in the literature so far. Some new stochastic properties of the process are derived and the corresponding optimal bivariate preventive maintenance policy with two decision parameters (age and operational history) is proposed. The latter is another novel feature of the study. Thus, development and application of the new mathematical models for modeling PM with the worse-than-minimal repair can be considered as the main contribution of the paper.

The structure of the paper is as follows: In Section 2, we introduce some preliminary results on the bivariate generalized Polya process (BVGPP) and the related repair process. In Section 3, we develop a bivariate preventive replacement policy assuming the BVGPP failure process and derive the corresponding long-run expected cost rate. In Section 4, we discuss the optimal policy for providing results of an illustrative numerical example. Finally, in Section 5, concluding remarks are given.

2. Preliminaries

In this section, we briefly review the definition of the bivariate generalized Polya process (BVGPP) and of some of its basic properties to be used in this paper. For this, we first need to recall the definition of the univariate generalized Polya process (GPP) via the concept of stochastic intensity. Note that, for an orderly (regular) counting process $\{N(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$ and its past history $H_{t-}\equiv \{N(u),\hspace{0.17em}\hspace{0.17em}0\le u<t\}$, the stochastic intensity is defined as (see, e.g., [13,22]),

$${\lambda }_t={\lim }_{\Delta t\to 0}\frac{P(N(t,t+\Delta t)=1|H_{t-})}{\Delta t}={\lim }_{\Delta t\to 0}\frac{E[N(t,t+\Delta t)|H_{t-}]}{\Delta t},$$

where $N(t_1,t_2)$, $t_1<t_2$, is the number of events in $[t_1,t_2)$. In the following definitions, $\lambda (t)$ is a non-negative deterministic function.

Definition 1 (Generalized Polya Process (GPP) [13]).

Let $\{N(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$be an orderly counting process and

(i)

$N(0)=0$;

(ii)

${\lambda }_t=(\alpha N(t-)+\beta )\lambda (t)$,

then it is called the Generalized Polya Process (GPP) with the corresponding parameter set$(\lambda (t),\alpha ,\beta )$,$\alpha \ge 0$,$\beta >0$.

As stated in [13], the GPP with $(\lambda (t),\alpha =0,\beta =1)$ reduces to the NHPP, and thus, the GPP is a generalization of the NHPP. Based on the GPP, and assuming that the repair times are negligible, Ref. [13] has defined a new type of imperfect repair for a system with the baseline (prior to the first repair) failure rate $\beta \lambda (t)$, which was called the ‘GPP repair’.

Definition 2 (GPP Repair).

If $\{N(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, where $N(t)$is the number of failures of the system in $(0,t]$, is the GPP with $(\lambda (t),\alpha ,\beta )$, then we say that the corresponding repair is the ‘GPP repair’ with the parameters $\alpha ,\beta >0$.

Accordingly, the corresponding stochastic intensity is given by

$${\lambda }_t=(\alpha N(t-)+\beta )\lambda (t)$$

(1)

Note that according to Definition 2 and Equation (1), the failure rate prior to the first failure starts from $\beta \lambda (t)$, which is called the baseline failure rate. From Equation (1), it is clear that the failure rate after each failure/repair is larger than that before it. Thus, due to GPP repair, the reliability performance of the system after failure/repair becomes worse. In general, in the definition of the GPP repair, the parameter $\beta$ can be set $\beta =1$ because the stochastic intensity in Equation (1) can be written as

$${\lambda }_t=(\frac{\alpha }{\beta }N(t-)+1)\beta \lambda (t)=({\alpha }^{\prime }N(t-)+1)\phi (t),$$

with ${\alpha }^{\prime }=\frac{\alpha }{\beta }$ and $\phi (t)=\beta \lambda (t)$. However, for a convenient description of the bivariate failure process, we follow Definition 2 throughout this paper.

Let $\{{\rm N}(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, where ${\rm N}(t)=(N_1(t),N_2(t))$, be a bivariate counting process and define the corresponding ‘pooled’ point process $\{M(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, where $M(t)=N_1(t)+N_2(t)$. The marginal point processes $\{N_i(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, for convenience, will be called type i point process, $i=1,2$, respectively. Furthermore, the events from type i point process $\{N_i(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$ will also be called type i events. For a regular multivariate process $\{{\rm N}(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, let $H_{Pt-}\equiv \{M(u),\hspace{0.17em}0\le u<t\}$ be the history of the pooled process in $[0,t)$, i.e., the set of all point events in $[0,t)$. Observe that $H_{Pt-}$ can equivalently be defined in terms of $M(t-)$ and the sequential arrival points of the events $0\le S_1\le S_2\le \cdots \le S_{M(t-)}<t$ in $[0,t)$, where $M(t-)$ is the total number of events in $[0,t)$ and $S_i$ is the time from 0 until the arrival of the ith event in $[0,t)$ of the pooled process $\{M(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$. Similarly, define the marginal histories of the marginal processes $H_{i,t-}\equiv \{N_i(u),\hspace{0.17em}0\le u<t\}$, $i=1,2$.

As with the case of univariate point processes, the most convenient general description of the multivariate point processes can be achieved through the stochastic intensities approach. Accord ingly, the ‘regular bivariate process’ can be specified by

$${\lambda }_{1t}\equiv \underset{\Delta t\to 0}{\lim }\frac{P(N_1(t,t+\Delta t)=1|H_{1,t-};H_{2,t-})}{\Delta t},$$

$${\lambda }_{2t}=\underset{\Delta t\to 0}{\lim }\frac{P(N_2(t,t+\Delta t)=1|H_{1,t-};H_{2,t-})}{\Delta t},$$

where $N_i(t_1,t_2),\hspace{0.17em}\hspace{0.17em}{t}_1<t_2$, denotes the number of events in $[t_1,t_2),\hspace{0.17em}i=1,2$, respectively (see [17]). According to [20], the BVGPP denoted further by BVGPP (${\lambda }_1(t),{\lambda }_2(t),\alpha ,\beta$), is defined as follows. In the following definition, ${\lambda }_i(t)$, $i=1,2$, are non-negative deterministic functions.

Definition 3 (Bivariate generalized Polya process (BVGPP)).

A bivariate counting process $\{{\rm N}(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$is called the bivariate generalized Polya process (BVGPP) with the set of parameters $({\lambda }_1(t),{\lambda }_2(t),\alpha ,\beta )$, ${\lambda }_i(t)\ge 0$for all $t\hspace{0.17em}\ge 0,\hspace{0.17em}\hspace{0.17em}i=1,2,\hspace{0.17em}\hspace{0.17em}\alpha \ge 0,\hspace{0.17em}\hspace{0.17em}\beta >0$, if

(i)

$N_1(0)=0,\hspace{0.17em}{N}_2(0)=0;$

(ii)

${\lambda }_{1t}=(\alpha (N_1(t-)+N_2(t-))+\beta ){\lambda }_1(t);$

(iii)

${\lambda }_{2t}=(\alpha (N_1(t-)+N_2(t-))+\beta ){\lambda }_2(t)$.

Conditions (ii) and (iii) in Definition 3 specify the dependence structure of the process in a fully intuitive way. That is, the occurrences of any type of events in the previous interval increase the occurrence probabilities of both types of events in the next interval. This type of dependency in a bivariate point process can be frequently observed in practice (see our examples in the Introduction).

Similar to the univariate counting process, a new type of dependent failure and repair process is defined based on the BVGPP in Definition 3, which is called ‘dependent worse-than-minimal repair process (DWMRP)’ [20]. Suppose that a system is composed of two parts (part 1 and part 2) having respective failure rates $\beta {\lambda }_i(t),\hspace{0.17em}i=1,2$. Under the DWMRP, the reliability performances of both parts after a repair of any part are worse than before the failure, which can be observed in reliability practice.

We will define now the concept of ‘thinning’ ([20,23]) for our further discussion and to provide some important properties of the BVGPP.

Definition 4 ($p(t)$-thinning).

Let$\{N(t),t\ge 0\}$be a univariate point process and denote it by $\{N_{p(\cdot )}(t),t\ge 0\}$the point process obtained by retaining (in the same location) every point of the process with probability $p(t)$and deleting it with probability $q(t)=1-p(t)$, independently of everything else. Denote by $\{N_{q(\cdot )}(t),t\ge 0\}$the point process constructed by the deleted points. Then the processes $\{N_{p(\cdot )}(t),t\ge 0\}$and $\{N_{q(\cdot )}(t),t\ge 0\}$are the $p(t)$-thinning of $\{N(t),t\ge 0\}$.

Denote:$\lambda (t)={\lambda }_1(t)+{\lambda }_2(t)$,$\Lambda (t)={\displaystyle {\int }_0^t\lambda (x)}dx$and$p_i(t)={\lambda }_i(t)/\lambda (t)$,$i=1,2$.

Proposition 1.

Let $\{N(t),t\ge 0\}$be the BVGPP$({\lambda }_1(t),{\lambda }_2(t),\alpha ,\beta )$. Then

(i)

$\{M(t),t\ge 0\}$is GPP$(\lambda (t),\alpha ,\beta )$.

(ii)

The process$\{{\rm N}(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$is constructed by$p_1(t)$-thinning of$\{M(t),t\ge 0\}$as$\{(M_{p_1(\cdot )}(t),\hspace{0.17em}{M}_{p_2(\cdot )}(t)),\hspace{0.17em}t\ge 0\}$.

(iii)

The marginal processes$\{N_i(t),t\ge 0\}$are GPP$({\gamma }_i(t),\alpha ,\beta )$, where

$${\gamma }_i(t)=\frac{{\lambda }_i(t)\exp \left\{\alpha {\displaystyle {\int }_0^t\hspace{0.17em}{\lambda }_1(x)+{\lambda }_2(x)dx}\right\}}{\alpha {\displaystyle {\int }_0^t{\lambda }_i(\nu )\exp \left\{\alpha {\displaystyle {\int }_0^{\nu }\hspace{0.17em}{\lambda }_1(x)+{\lambda }_2(x)dx}\right\}d\nu }+1},\hspace{0.17em}\hspace{0.17em}i=1,2.$$

See [20] for the proof of Proposition 1. The following proposition presents the joint distribution of number of events ([20]).

Proposition 2.

Let $t>0$. It holds that

$$\begin{array}{l}P(N_i(t)=n_i,i=1,2)\\ =\frac{\Gamma (\beta /\alpha +n_1+n_2)}{\Gamma (\beta /\alpha )n_1!n_2!}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{n_1}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{n_2}\\ \times {(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }.\end{array}$$

3. Bivariate Preventive Replacement Policy

We will now develop a new preventive replacement policy for a repairable deteriorating system which is composed of two statistically dependent parts. It should be noted that the PM models based on the univariate counting processes were only considered in the literature previously. Denote by $N_i(t)$, $i=1,2$, the number of failures in part 1 and part 2 until time $t$, respectively. Under the BVGPP failure/repair process, we assume that $\{{\rm N}(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$, where ${\rm N}(t)=(N_1(t),N_2(t))$, follows BVGPP (${\lambda }_1(t),{\lambda }_2(t),\alpha ,\beta$). As mentioned before, under the BVGPP (or DWMRP), the reliability performances of ‘both parts’ deteriorate on each failure of any of the two parts, as the corresponding stochastic intensities in Definition 3 ‘count’ the overall number of events, i.e., $N_1(t-)+N_2(t-)$. Therefore, it could be reasonable to suggest the preventive replacement policy based on $N_1(t)+N_2(t)$ for the BVGPP $({\lambda }_1(t),{\lambda }_2(t),\alpha ,\beta )$ failure/repair process. Recall that $S_N$ ($N=1,2,\dots$) denotes the arrival times in the pooled point process $\{M(t),\hspace{0.17em}\hspace{0.17em}t\ge 0\}$.

3.1. Bivariate Preventive Replacement Policy

The system is replaced at time $T$$(T>0)$ or at $S_N$ ($N=1,2,\dots$) after its inception into operation (or last replacement), whichever occurs first, and it undergoes the BVGPP repairs at failures between replacements. The times for repairs and replacements are negligible.

Let us denote, by $c(N,T)$, the corresponding long-run expected cost rate function. Let $c_{GPP}^{(i)}$ be the cost incurred by a BVGPP repair performed on the failure of part i, $i=1,2$, and $c_r$ be the cost of system’s replacement. To derive the cost rate function, we need some preliminary lemmas. In the following, denote by $N_{\hspace{0.17em}GP}^{\hspace{0.17em}\hspace{0.17em}(i)},\hspace{0.17em}\hspace{0.17em}i=1,2$, the total number of BVGPP repairs of part $i$ in a renewal cycle (between replacements).

Lemma 1.

Conditional expectations $E(N_{GP}^{(i)}|S_N\le T)$, i = 1,2, are given by

$$\begin{array}{ll}E(N_{GP}^{(i)}|S_N\le T)& ={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_i(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_i(t)}{\lambda (t)}\right\}}\\ & \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\cdot \frac{1}{P(S_N\le T)},i=1,2. \end{array}$$

Proof. 

In this proof, we derive just $E(N_{GP}^{(1)}|S_N\le T)$, whereas $E(N_{GP}^{(2)}|S_N\le T)$ can be obtained ‘symmetrically’. Observe that

$$E(N_{GP}^{(1)}|S_N\le T)={\displaystyle {\int }_0^{T}E(N_{GP}^{(1)}|S_N\le T,S_N=t)f_{(S_N|S_N\le T)}(t)}dt,$$

(2)

where $f_{(S_N|S_N\le T)}(t)$ is the conditional pdf of $(S_N|S_N\le T)$ given by

$$\begin{array}{ll}{f}_{(S_N|S_N\le T)}(t)& =\frac{f_{S_N}(t)}{P(S_N\le T)}\\ & =\frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)\cdot \frac{1}{P(S_N\le T)}. \end{array}$$

Furthermore,

$$\begin{array}{ll}E(N_{GP}^{(1)}|S_N\le T,S_N=t)& =E(N_{GP}^{(1)}|S_N=t,I_N=1)P(I_N=1|S_N=t)\\ & +E(N_{GP}^{(1)}|S_N=t,I_N=2)P(I_N=2|S_N=t),\end{array}$$

(3)

where $I_N=i$, $i=1,2$, if the failure at time $S_N$ occurs in part $i$, respectively. From Proposition 1-(ii),

$$P(I_N=i|S_N=t)={\lambda }_i(t)/\lambda (t),\hspace{0.17em}\hspace{0.17em}i=1,2,$$

(4)

and $E(N_{GP}^{(1)}|S_N=t,I_N=1)$ can be represented as

$$E(N_{GP}^{(1)}|S_N=t,I_N=1)={\displaystyle \sum _{n_1=1}^N{n}_{1}P(N_{GP}^{(1)}=n_1|S_N=t,I_N=1)}={\displaystyle \sum _{n_1=1}^N{n}_1\frac{f_{(N_{GP}^{(1)},S_N,I_N)}(n_1,t,1)}{f_{(S_N,I_N)}(t,1)}}$$

where $f_{(S_N,I_N)}(t,1)$ is the joint distribution of $(S_N=t,I_N=1)$ and $f_{(N_{GP}^{(1)},S_N,I_N)}(n_1,t,1)$ is that of $(N_{GP}^{(1)}=n_1,S_N=t,I_N=1)$, given by

$$\begin{array}{l}{f}_{(S_N,I_N)}(t,1)\\ =\frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)\frac{{\lambda }_1(t)}{\lambda (t)}, \end{array}$$

and

$$\begin{array}{l}{f}_{(N_{GP}^{(1)},S_N,I_N)}(n_1,t,1)\\ ={\lim }_{\Delta t\to 0}\frac{1}{\Delta t}P(N_1(t-)=n_1-1,\hspace{0.17em}{N}_2(t-)=N-n_1,t\le S_N<t+\Delta t,I_N=1)\\ =P(N_1(t-)=n_1-1,\hspace{0.17em}{N}_2(t-)=N-n_1)\\ \times {\lim }_{\Delta t\to 0}\frac{1}{\Delta t}P(t\le S_N<t+\Delta t,I_N=1|N_1(t-)=n_1-1,\hspace{0.17em}{N}_2(t-)=N-n_1)\\ =\frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(n_1-1)!(N-n_1)!}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{n_1-1}\\ \times {\left(\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{N-n_1}\exp {\{-\alpha \Lambda (t)\}}^{\beta /\alpha }\\ \times ((N-1)\alpha +\beta )\lambda (t)\frac{{\lambda }_1(t)}{\lambda (t)}.\end{array}$$

Therefore, we have

$$\begin{array}{l}P(N_{GP}^{(1)}=n_1|S_N=t,I_N=1)\\ =\frac{(N-1)!}{(n_1-1)!(N-n_1)!}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{n_1-1}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{N-n_1}\\ \raisebox{1ex}{$$}\!\left/ \!\raisebox{-1ex}{$$}\right.{(1-\exp \{-\alpha \Lambda (t)\})}^{N-1}\\ =\frac{(N-1)!}{(n_1-1)!(N-n_1)!}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}}{1-\exp \{-\alpha \Lambda (t)\}}\right)}^{n_1-1}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}}{1-\exp \{-\alpha \Lambda (t)\}}\right)}^{N-n_1}\\ =\frac{(N-1)!}{(n_1-1)!(N-n_1)!}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^{n_1-1}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^{N-n_1}, n_1=1,2,\dots ,N. \end{array}$$

Let $L\equiv N_{GP}^{(1)}-1$. Then, the conditional distribution of $P(L=m|S_N=t,I_N=1)$, $m=0,1,2,\dots ,N-1$, is

$$\begin{array}{l}P(L=m|S_N=t,I_N=1)=P(N_{GP}^{(1)}=m+1|S_N=t,I_N=1)\\ =\frac{(N-1)!}{m!(N-1-m)!}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^m{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^{N-1-m},\end{array}$$

$m=0,1,2,\dots ,N-1$, which is the Binomial distribution with parameters $N-1$ and $\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}$. Accordingly, $E(N_{GP}^{(1)}|S_N=t,I_N=1)$ is given by

$$\begin{array}{l}E(N_{GP}^{(1)}|S_N=t,I_N=1)=E(L+1|S_N=t,I_N=1)\\ =(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+1. \end{array}$$

In a similar way,

$$\begin{array}{l}{f}_{(S_N,I_N)}(t,2)\\ =\frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)\frac{{\lambda }_2(t)}{\lambda (t)}, \end{array}$$

and

$$\begin{array}{l}{f}_{(N_{GP}^{(1)},S_N,I_N)}(n_1,t,2)\\ ={\lim }_{\Delta t\to 0}\frac{1}{\Delta t}P(N_1(t-)=n_1,\hspace{0.17em}{N}_2(t-)=N-1-n_1,t\le S_N<t+\Delta t,I_N=2)\\ =P(N_1(t-)=n_1,\hspace{0.17em}{N}_2(t-)=N-1-n_1)\\ \times {\lim }_{\Delta t\to 0}\frac{1}{\Delta t}P(t\le S_N<t+\Delta t,I_N=2|N_1(t-)=n_1,\hspace{0.17em}{N}_2(t-)=N-1-n_1)\\ =\frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )n_1!(N-1-n_1)!}{\left(\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{n_1}\\ \times {\left(\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{-\alpha [\Lambda (t)-\Lambda (x)]\}dx}\right)}^{N-1-n_1}\exp {\{-\alpha \Lambda (t)\}}^{\beta /\alpha }\\ \times ((N-1)\alpha +\beta )\lambda (t)\frac{{\lambda }_2(t)}{\lambda (t)}. \end{array}$$

Thus,

$$\begin{array}{l}P(N_{GP}^{(1)}=n_1|S_N=t,I_N=2)\\ =\frac{(N-1)!}{n_1!(N-1-n_1)!}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^{n_1}{\left(\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right)}^{N-1-n_1},\end{array}$$

$n_1=0,1,\dots ,N-1$, and

$$E(N_{GP}^{(1)}|S_N=t,I_N=2)=(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}.$$

From Equations (3) and (4),

$$\begin{array}{ll}E(N_{GP}^{(1)}|S_N\le T,S_N=t)& =\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}\right\}\left(\frac{{\lambda }_1(t)}{\lambda (t)}+\frac{{\lambda }_2(t)}{\lambda (t)}\right)+\frac{{\lambda }_1(t)}{\lambda (t)}\\ & =(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_1(t)}{\lambda (t)}, \end{array}$$

and from Equation (2),

$$\begin{array}{l}E(N_{GP}^{(1)}|S_N\le T)={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_1(t)}{\lambda (t)}\right\}}\\ \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\cdot \frac{1}{P(S_N\le T)}. \end{array}$$

□

Lemma 2.

$E(N_{GP}^{(i)})$and $i=1,2$, aregiven by

$$\begin{array}{l}E(N_{GP}^{(1)})={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_1(t)}{\lambda (t)}\right\}}\\ \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\\ +{\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1-j}j\cdot \frac{\Gamma (\beta /\alpha +j+k)}{\Gamma (\beta /\alpha )j!k!}{\left(\alpha {\displaystyle {\int }_0^T{\lambda }_1(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^j}}\\ \times {\left(\alpha {\displaystyle {\int }_0^T{\lambda }_2(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^k{(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }, \end{array}$$

(5)

and

$$\begin{array}{l}E(N_{GP}^{(2)})={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_2(t)}{\lambda (t)}\right\}}\\ \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\\ +{\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1-j}j\cdot \frac{\Gamma (\beta /\alpha +j+k)}{\Gamma (\beta /\alpha )j!k!}{\left(\alpha {\displaystyle {\int }_0^T{\lambda }_2(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^j}}\\ \times {\left(\alpha {\displaystyle {\int }_0^T{\lambda }_1(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^k{(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }.\end{array}$$

(6)

Proof. 

Observe that, using Proposition 2,

$$\begin{array}{l}E(N_{GP}^{(1)}|S_N>T)=E(N_1(T)|M(T)\le N-1)\\ ={\displaystyle \sum _{j=0}^{N-1}j\hspace{0.17em}P(N_1(T)=j|M(T)\le N-1)}={\displaystyle \sum _{j=0}^{N-1}j\hspace{0.17em}\frac{P(N_1(T)=j,N_2(T)\le N-1-j\hspace{0.17em})}{P(M(T)\le N-1)}}\\ ={\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1-j}j\cdot \frac{\Gamma (\beta /\alpha +j+k)}{\Gamma (\beta /\alpha )j!k!}{\left(\alpha {\displaystyle {\int }_0^T{\lambda }_1(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^j}}\\ \times {\left(\alpha {\displaystyle {\int }_0^T{\lambda }_2(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^k{(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }\cdot \frac{1}{P(S_N>T)}. \end{array}$$

Therefore, using the result in Lemma 1, we have

$$\begin{array}{ll}E(N_{GP}^{(1)})& =E(N_{GP}^{(1)}|S_N\le T)P(S_N\le T)+E(N_{GP}^{(1)}|S_N>T)P(S_N>T)\\ & ={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_1(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_1(t)}{\lambda (t)}\right\}}\\ & \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\\ & +{\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1-j}j\cdot \frac{\Gamma (\beta /\alpha +j+k)}{\Gamma (\beta /\alpha )j!k!}{\left(\alpha {\displaystyle {\int }_0^T{\lambda }_1(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^j}}\\ & \times {\left(\alpha {\displaystyle {\int }_0^T{\lambda }_2(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^k{(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }, \end{array}$$

and symmetrically, $E(N_{GP}^{(2)})$ can be derived as

$$\begin{array}{ll}E(N_{GP}^{(2)})& =E(N_{GP}^{(2)}|S_N\le T)P(S_N\le T)+E(N_{GP}^{(2)}|S_N>T)P(S_N>T)\\ & ={\displaystyle \underset{0}{\overset{T}{\int }}\left\{(N-1)\frac{\alpha {\displaystyle {\int }_0^t{\lambda }_2(x)\exp \{\alpha \Lambda (x)\}dx}}{\exp \{\alpha \Lambda (t)\}-1}+\frac{{\lambda }_2(t)}{\lambda (t)}\right\}}\\ & \times \frac{\Gamma (\beta /\alpha +N-1)}{\Gamma (\beta /\alpha )(N-1)!}{\left(1-\exp \{-\alpha \Lambda (t)\}\right)}^{N-1}\exp \{-\beta \Lambda (t)\}((N-1)\alpha +\beta )\lambda (t)dt\\ & +{\displaystyle \sum _{j=0}^{N-1}{\displaystyle \sum _{k=0}^{N-1-j}j\cdot \frac{\Gamma (\beta /\alpha +j+k)}{\Gamma (\beta /\alpha )j!k!}{\left(\alpha {\displaystyle {\int }_0^T{\lambda }_2(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^j}}\\ & \times {\left(\alpha {\displaystyle {\int }_0^T{\lambda }_1(x)\exp \{-\alpha [\Lambda (T)-\Lambda (x)]\}dx}\right)}^k{(\exp \{-\alpha \Lambda (T)\})}^{\beta /\alpha }. \end{array}$$

□

In the following theorem, we derive the corresponding expected cost rate function $c(N,T)$, which is, as usual, defined as the expected cost on a renewal cycle over the expected length of this cycle.

Theorem 1.

The cost rate function $c(N,T)$is given by

$$\begin{gathered} c(N,T)=\frac{\left(c_{GPP}^{(1)}E[N_{GP}^{(1)}]+c_{GPP}^{(2)}E[N_{GP}^{(2)}]+c_r\right)}{{\displaystyle \underset{0}{\overset{T}{\int }}\hspace{0.17em}{\displaystyle \sum _{j=0}^{N-1}\frac{\Gamma (\beta /\alpha +j)}{\Gamma (\beta /\alpha )j!}{\left(1-\exp \left\{-\alpha \Lambda (t)\right\}\right)}^j{\left(\exp \left\{-\alpha \Lambda (t)\right\}\right)}^{\frac{\beta }{\alpha }}}}\hspace{0.17em}dt}, \\ N=1,2,\dots , T>0, \end{gathered}$$

(7)

where$E[N_{GP}^{(i)}]$,$i=1,2$, are given by Equations (5) and (6), respectively.

Proof. 

Observe that the expected length of one renewal cycle is $E[\min (S_N,T)]$. Then, by the Renewal Reward Theorem ([23]) the long-run expected cost rate function $c(N,T)$ is defined as

$$c(N,T)=\frac{c_{GPP}^{(1)}E[N_{GP}^{(1)}]+c_{GPP}^{(2)}E[N_{GP}^{(2)}]+c_r}{E[\min (S_N,T)]}.$$

(8)

From Proposition 1-(i), the expected length of a cycle, $E[\min (S_N,T)]$, is given as follows:

$$\begin{array}{ll}E[\min (S_N,T)]& ={\displaystyle \underset{0}{\overset{T}{\int }}P(\min (S_N,T)>t)}\hspace{0.17em}dt={\displaystyle \underset{0}{\overset{T}{\int }}P(S_N>t)}\hspace{0.17em}dt={\displaystyle \underset{0}{\overset{T}{\int }}P(M(t)\le N-1)}\hspace{0.17em}dt\\ & ={\displaystyle \underset{0}{\overset{T}{\int }}\hspace{0.17em}{\displaystyle \sum _{j=0}^{N-1}\frac{\Gamma (\beta /\alpha +j)}{\Gamma (\beta /\alpha )j!}{\left(1-\exp \left\{-\alpha \Lambda (t)\right\}\right)}^j{\left(\exp \left\{-\alpha \Lambda (t)\right\}\right)}^{\frac{\beta }{\alpha }}}}\hspace{0.17em}dt.\end{array}$$

(9)

Eventually, using Lemma 2, Equations (8) and (9), the long-run expected cost rate function $c(N,T)$ is given by Equation (7). □

It is important to compare the proposed bivariate policy (age or the occurrence of the N-th event in the corresponding BVGGP, whichever comes first) with the conventional ordinary age-based replacement policy at age $T_A\hspace{0.17em}(T_A>0)$. Thus, the number of GPP repairs are not considered as an additional parameter in this simplified policy. The corresponding expected cost rate (denoted by $c_A(T_A)\hspace{0.17em}$) can be obtained as

$$c_A(T_A\hspace{0.17em})=\frac{\frac{\beta }{\alpha }\left\{c_{GPP}^{(1)}{\displaystyle {\int }_0^{T_A}\alpha {\lambda }_1(x)\exp (\alpha \Lambda (x))dx}+c_{GPP}^{(2)}{\displaystyle {\int }_0^{T_A}\alpha {\lambda }_2(x)\exp (\alpha \Lambda (x))dx}\right\}+c_r}{T_A}.$$

(10)

Observe that Equation (10) can be directly derived from Theorem 1 by setting $N=\infty$, i.e., $c_A(T_A)\hspace{0.17em}=\underset{N\to \infty }{\lim }c(N,T_A)$. In addition, we can see that when $c_{GPP}^{(1)}=c_{GPP}^{(2)}=c_{GPP}$, the expected cost rate $c_A(T_A)\hspace{0.17em}$ in Equation (10) is equal to that in [13]:

$$c{*}_A(T_A\hspace{0.17em})=\frac{c_{GPP}\beta \left\{\exp (\alpha \Lambda (T_A))-1\right\}/\alpha +c_r}{T_A}$$

3.2. Optimal PM

We can now formulate the optimal PM problem for the described setting. Thus, the optimal vector $(N^{*},T^{*})$ should be obtained such that

$$c(N^{*},T^{*})={\min }_{T>0,N=1,2,\dots }c(N,T).$$

To find $(N^{*},T^{*})$, the two-stage procedure will be applied. At the first step, for a fixed $N$, we find $T^{*}(N)$ such that

$$c(N,T^{*}(N))={\min }_{T>0}c(N,T).$$

At the second step, we search for $N^{*}$ such that

$$c(N^{*},T^{*}(N^{*}))={\min }_{N=1,2,\dots }c(N,T^{*}(N)).$$

Then, the optimal maintenance policy parameters are given by $(N^{*},T^{*}(N^{*}))$.

The expression for $c(N,T)$ obtained in Theorem 1 is extremely cumbersome and its analytical analysis of the optimal solution is practically impossible. Therefore, in the next section, we will illustrate our findings numerically. Note that, from general considerations, it is clear that the optimal policy proposed in this study should result in a smaller (not larger) optimal expected cost rate than for the case defined by Equation (10). The numerical study of the next section among other findings illustrates this claim as well.

4. Numerical Illustration and Discussion

Following the optimization procedure stated above, we conduct numerical studies for illustration. Suppose that two parts of the system has the following baseline intensities: ${\lambda }_1(t)=0.25\hspace{0.17em}(t+2)$ and ${\lambda }_2(t)=0.5\hspace{0.17em}(t+2)$. Let the repair and replacement costs are given by $c_{GPP}^{(1)}=5$, $c_{GPP}^{(2)}=10$, and $c_r=100$, respectively. Then, for instance, for $\alpha =\beta =1$, the optimal values are obtained as $(N^{*},T^{*})=(5,\hspace{0.17em}3.9801)$ and the corresponding minimal cost rate is $c(N^{*},T^{*})=130.8159$.

In what follows, assume $\beta =1$, which is a natural assumption in defining the BVGGP that describes the failure/repair process. Then the optimal preventive maintenance policy $(N^{*},T^{*})$ and the corresponding cost rate $c(N^{*},T^{*})$ for different values $\hspace{0.17em}{c}_{GPP}^{(1)},\hspace{0.17em}{c}_{GPP}^{(2)},$ and $\alpha$ are given in

Table 1. The optimal $(N^{*},T^{*})$ and $c(N^{*},T^{*})$ for different values of $\alpha$ and $c_{GPP}^{(i)},i=1,2$ ($c_r=80$).

	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(1)}$	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(2)}$	${\mathit{N}}^{*}$	${\mathit{T}}^{*}$	$\mathit{E}[\mathbf{min}({\mathit{S}}_{{\mathit{N}}^{*}},{\mathit{T}}^{*})]$	$\mathit{c}({\mathit{N}}^{*},{\mathit{T}}^{*})$
$\alpha =0.1$	1	1	43	5.167063	4.839603	24.613510
		5	25	3.570791	3.515909	40.370606
		10	22	2.917147	2.907530	54.556764
	5	1	29	4.090193	3.966627	33.444139
		5	23	3.245566	3.217194	46.421183
		10	22	2.757668	2.753356	59.567110
	10	1	24	3.481198	3.429802	41.947039
		5	22	2.963350	2.951428	53.259190
		10	22	2.598863	2.597111	65.502361
$\alpha =0.5$	1	1	27	3.813358	2.381693	44.905501
		5	12	2.643918	1.847879	66.290938
		10	10	1.944151	1.622081	84.952145
	5	1	16	2.983852	2.046325	57.014955
		5	12	2.148094	1.77202	74.384448
		10	10	1.771289	1.556203	91.553627
	10	1	12	2.480472	1.832928	68.366798
		5	10	1.998873	1.638653	83.243218
		10	9	1.683485	1.475049	99.086616
$\alpha =1.0$	1	1	23	4.458552	1.694946	60.768793
		5	10	3.085826	1.365100	85.423285
		10	7	2.371505	1.204417	106.546497
	5	1	13	3.641256	1.473804	74.858730
		5	8	2.904525	1.268782	94.502131
		10	6	2.364720	1.134614	114.03389
	10	1	9	3.205668	1.320562	87.817809
		5	7	2.499525	1.206568	104.632805
		10	6	1.960038	1.120323	122.666147

and

Table 2. The optimal $(N^{*},T^{*})$ and $c(N^{*},T^{*})$ for different values of $\alpha$ and $c_{GPP}^{(i)},i=1,2$ ($c_r=100$).

	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(1)}$	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(2)}$	${\mathit{N}}^{*}$	${\mathit{T}}^{*}$	$\mathit{E}[\mathbf{min}({\mathit{S}}_{{\mathit{N}}^{*}},{\mathit{T}}^{*})]$	$\mathit{c}({\mathit{N}}^{*},{\mathit{T}}^{*})$
$\alpha =0.1$	1	1	48	5.485968	5.062683	28.648590
		5	27	3.817166	3.737695	45.883026
		10	23	3.129246	3.111605	61.198240
	5	1	32	4.357206	4.196250	38.343312
		5	24	3.481181	3.429788	52.433799
		10	22	2.963349	2.951427	66.573986
	10	1	26	3.722206	3.649914	47.592342
		5	23	3.178678	3.156975	59.803276
		10	22	2.794272	2.789043	72.923917
$\alpha =0.5$	1	1	31	4.078488	2.464836	53.138310
		5	14	2.779690	1.953065	76.74202
		10	12	2.014106	1.730446	97.289067
	5	1	18	3.240933	2.127228	66.583640
		5	12	2.480472	1.832928	85.458497
		10	10	1.998873	1.638653	104.054023
	10	1	13	2.757591	1.905457	79.035012
		5	12	2.069472	1.749067	95.393134
		10	10	1.808951	1.572418	112.421272
$\alpha =1.0$	1	1	27	4.639853	1.754148	72.399750
		5	11	3.501533	1.405433	99.841511
		10	8	2.563363	1.266254	122.837466
	5	1	15	3.878627	1.530984	88.176667
		5	9	3.205668	1.320562	109.772261
		10	7	2.499525	1.206568	130.791007
	10	1	11	3.176349	1.405078	102.448963
		5	8	2.704431	1.267669	120.746421
		10	7	2.06428	1.193740	140.30300

for $c_r=80$ and $c_r=100$, respectively. From these tables, we can observe that as the degree of the GPP repair increases (i.e., $\alpha$ increases) and as each GPP repair cost incurred by the failure of each part increases (i.e., $c_{GPP}^{(i)}$ increases), the mean time until replacement, $E[\min (S_{N^{*}},T^{*})]$, decreases; that is, the system should be replaced earlier. Moreover, it can be seen that as the replacement cost $c_r$ gets larger, the mean time until replacement, $E[\min (S_{N^{*}},T^{*})]$, increases.

To compare the cost rates of the two policies, $c(N^{*},T^{*})$ and $c_A(T_A{}^{*})$ (see our discussion at the end of Section 3), we introduce the following index:

$$\Delta (\%)=\frac{c_A(T_A^{*})-c(N^{*},T^{*})}{c_A(T_A^{*})}\times 100\hspace{0.17em}\%$$

which indicates the relative difference between the minimum cost rate of the two maintenance policies. A larger $\Delta \hspace{0.17em}(\%)$ means that the proposed policy has priority over the conventional age-based replacement policy from the cost-rate-minimization point of view. The optimal conventional maintenance policy $T_A{}^{*}$ and the corresponding cost rates $c_A(T_A{}^{*})$ under different combinations of $c_r,\hspace{0.17em}\hspace{0.17em}{c}_{GPP}^{(1)},\hspace{0.17em}{c}_{GPP}^{(2)},\hspace{0.17em}$ and $\alpha$ are summarized in

Table 3. The optimal $T_A{}^{*}$, $c(T_A{}^{*})$, and $\Delta (\%)$ for different values of $\alpha$ and $c_{GPP}^{(i)},i=1,2$ ($c_r=80$ and $c_r=100$).

	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(1)}$	${\mathit{c}}_{\mathit{G}\mathit{P}\mathit{P}}^{(2)}$		${\mathit{c}}_{\mathit{r}}=80$			${\mathit{c}}_{\mathit{r}}=100$
			${\mathit{T}}_{\mathit{A}}{}^{*}$	$\mathit{c}({\mathit{T}}_{\mathit{A}}{}^{*})$	$\mathbf{\Delta }(\mathit{\%})$	${\mathit{T}}_{\mathit{A}}{}^{*}$	$\mathit{c}({\mathit{T}}_{\mathit{A}}{}^{*})$	$\mathbf{\Delta }(\mathit{\%})$
$\alpha =0.1$	1	1	4.792556	24.73827	0.504320	5.022464	28.81093	0.563467
		5	3.497507	40.41667	0.113973	3.712159	45.96080	0.169218
		10	2.902790	54.56740	0.019491	3.102907	61.22304	0.040508
	5	1	3.936392	33.52346	0.236613	4.158405	38.46142	0.307082
		5	3.206681	46.44719	0.055993	3.414981	52.48364	0.094965
		10	2.750663	59.57278	0.009518	2.946040	66.58946	0.023238
	10	1	3.414981	41.98691	0.094961	3.627975	47.66229	0.146758
		5	2.946040	53.27157	0.023239	3.147431	59.83133	0.046888
		10	2.595846	65.50489	0.003861	2.785933	72.93223	0.011398
$\alpha =0.5$	1	1	2.280309	47.06474	4.587806	2.370803	55.65973	4.530062
		5	1.75982	69.16807	4.159625	1.847788	80.2483	4.304914
		10	1.510966	87.83091	3.277622	1.595671	100.6977	3.385016
	5	1	1.938804	59.70407	4.470884	2.028205	69.78078	4.568952
		5	1.639218	77.22535	3.678717	1.725783	89.10439	4.09171
		10	1.445795	94.25774	2.86885	1.529378	107.693	3.379028
	10	1	1.725783	71.28351	4.091706	1.813389	82.57799	4.290463
		5	1.529378	86.15439	3.379018	1.614373	98.86931	3.51593
		10	1.378737	101.7853	2.651349	1.460969	115.8614	2.969175
$\alpha =1.0$	1	1	1.550914	67.21102	9.585075	1.609079	79.86217	9.344124
		5	1.215585	95.30571	10.36918	1.272421	111.3726	10.35361
		10	1.05409	118.2154	9.870882	1.109187	136.6934	10.13651
	5	1	1.33114	83.43543	10.27925	1.388719	98.13309	10.14577
		5	1.137439	105.2673	10.22651	1.193554	122.416	10.3285
		10	1.011592	125.9852	9.486281	1.06608	145.2245	9.938745
	10	1	1.193554	97.93283	10.32853	1.250203	114.2905	10.36091
		5	1.06608	116.1796	9.938746	1.121332	134.4538	10.19486
		10	0.967724	135.0222	9.151127	1.021496	155.117	9.55021

. From Table 1, Table 2 and Table 3, we can see that under all combinations $c_r,\hspace{0.17em}\hspace{0.17em}{c}_{GPP}^{(1)},\hspace{0.17em}{c}_{GPP}^{(2)},$ and $\alpha$, the cost rate for the proposed policy is relatively smaller than that for the original one having the values of $\Delta \hspace{0.17em}(\%)$ in a range of (10%, 20%). This means that there exists a meaningful difference between the two maintenance policies in terms of the expected cost rate. Table 3 shows that the difference between the two minimum cost rates decreases as $\alpha$ increases.

The graphs for the optimal $N^{*}$ and $T^{*}$ with respect to the value of $\alpha$ are given in Figure 1 and Figure 2. As $c_{GPP}^{(1)}$ or $c_{GPP}^{(2)}$ increases, the system should be replaced earlier and thus the corresponding curves are ordered. As $\alpha$ increases, the mean length of the renewal cycle $E[\min (S_{N^{*}},T^{*})]$ should be smaller and thus $N^{*}$ and $T^{*}$ initially decreases. However, when $\alpha$ is larger, the replacement should be made mainly based on the number of failures $N^{*}$ (this follows from the form of stochastic intensities in Definition 3) and the role of $T^{*}$ should be weaker. Due to this effect, $T^{*}$ is increasing when $\alpha$ is increasing.

In Figure 3, the minimum expected cost rates of the two policies for different values of $\alpha$ are presented. The superiority of the proposed replacement policy is clearly seen.

5. Concluding Remarks

In reliability modeling and analysis, the assumption of independence of components in a system is usually made for simplicity and convenience of stochastic description. However, most often, the failures in two or more parts in a system are statistically dependent. Furthermore, in practice, a failure of one part in a system often causes additional stress or damage to the remaining parts, which results in a worse condition of a system than it had just prior this failure. Even a minimal repair of the failed component results in a worse-than-minimal repair of a system in this case.

To model the described practical setting, we employ the bivariate generalized Polya process, which corresponds to the dependent worse-than-minimal repair process. Under these assumptions, a new bivariate preventive maintenance policy based on two parameters (age and operational history) has been proposed and discussed. The corresponding long-run average cost rate has been derived and the optimal replacement policies are investigated and illustrated numerically.

Our mathematical study has a clear practical application in the field of the PM modeling. Along with optimal maintenance, in the future research, the proposed concept of dependent worse-than-minimal repair processes could be used also for describing other reliability properties of repairable systems, such as stationary and non-stationary availability.

As far as we know, our study is the first to apply the dependent bivariate or multivariate counting processes to modeling the multivariate failure processes for stochastic description of repairable systems. Some new stochastic properties of these processes have been derived and the corresponding optimal bivariate preventive maintenance policy with two decision parameters (age and operational history) has been proposed. The latter is another novel feature of the study. Thus, development and application of the new mathematical models for modeling PM with the worse-than-minimal repair can be considered as the main contribution of the paper.

The developed approach and obtained results provide the tools for more adequate stochastic description of real systems with dependent components and worse-than-minimal repair. Neglecting these real-life properties can result in substantial discrepancies in reliability estimates and PM schedules of systems, along with higher costs of maintenance.

The implications of some assumptions in the study should be also addressed in future research. In the multivariate setting, it is also interesting to consider the case when some components undergo minimal repair, whereas others undergo worse-than-minimal repair.