# Condition-based maintenance for a system subject to multiple degradation
processes with stochastic arrival intensity

Lucía Bautista  Department of Mathematics
University of Extremadura, Cáceres (Spain)
email<EMAIL_ADDRESS>Inma T. Castro  Department of Mathematics
University of Extremadura, Cáceres (Spain)
email<EMAIL_ADDRESS>Luis Landesa  Department of Computers and
Communications Technology
University of Extremadura, Cáceres (Spain)
email<EMAIL_ADDRESS>
###### Abstract

In this work, a system subject to different deterioration processes is
analysed. The arrival of the degradation processes to the system is modelled
using a shot-noise Cox process. The degradation processes grow according to an
homogeneous gamma process. The system fails when a degradation process exceeds
a failure threshold. The combined process of initiation and growth of the
degradation processes is modelled and the system reliability is obtained.
Heterogeneities are  also integrated in the model assuming that the  inverse
of the scale parameter follows a uniform distribution.  A maintenance strategy
is implemented in this system and the state of the system is checked in
inspection times. If the system is  working at inspection time, a preventive
replacement is performed if the deterioration level of a degradation process
exceeds a certain threshold. A corrective replacement is performed if the
system  is down at inspection time. Under this maintenance strategy, the
expected cost rate is obtained.  Sensitivity analysis on the main parameters
of the gamma process is performed.

###### keywords:

condition-based maintenance, gamma process, heterogeneities, shot-noise Cox
process, complex degradation processes.

##  1 Introduction

It is well known that Condition-Based Maintenance (CBM) is a maintenance
strategy based on monitoring the operating condition of a system. Compared to
time-based maintenance, and thanks to the development of sensor technologies,
CBM usually results in lower maintenance costs, avoiding unnecessary
preventive maintenance activities and reducing unexpected failures [21].
However, as Alaswad and Xiang [1] claim, further research on CBM is still in
great need. Since industrial systems are becoming  increasingly complex and
they are likely to suffer from multiple degradation processes, further
research is needed in modelling the maintenance when a system possesses more
than one degradation path. Systems subject to multiple degradation processes
can be found in electronic products, heavy machine tools and piping systems
[28]. On a pavement network, several different degradation processes, such as
fatigue cracking and pavement deformation, may develop simultaneously [39].

In literature, significant and meaningful prior research has been done on
reliability and maintenance policies for systems subject to multiple
degradation processes. An earlier work on this topic is [16], where a multiple
multivariate regression is applied to model the complexity of the processes.
Some studies model the joint probability assuming that the degradation
processes are independent ([6], [5], [4]). It could be the case of components
in different production units without interference. But, in some cases,
interaction between adjacent processes has a significant influence on the
propagation characteristics of the rest of the processes. In these cases, the
dependence between degradation processes has to be taken into account.
Different works that deal with dependent degradation processes can be found in
[38] and [8].

In general, for complex degradation processes, there are also two approaches
to model the degradation mechanism of the processes. The first approach
considers that all the processes start to degrade at the same time ([23]).
However, as Kuniewski et. al claim [24], it is unlikely that all processes
appear at the same time. The second approach assumes that the processes
initiate at random times and then grow depending on the environment and
conditions of the system. In this approach, two stochastic processes have to
be combined: the initiation process and the growth process. Let
$\left\\{N^{*}(t),t\geq 0\right\\}$ be the point process of the degradation
processes with arrival times $S_{1}<S_{2}<\ldots$ (“initiating events”). Each
degradation process triggers the “effective event” of failure after a random
time $D_{i}$, where $D_{i}$ are i.i.d. non-negative variables. The sequence
$\left\\{S_{i}+D_{i}\right\\}$ forms a new point process. The initiating
events can often be interpreted as potentially harmful events affecting a
system. The “wear process” triggered is activated at the moment of the
occurrence of the harmful event and continuously increases with time [9]. This
approach has been recently dealt by Cha and Filkenstein [11] assuming that the
initiation process follows a generalized Polya process. They obtained the
survival function of the system and an analysis of the residual lifetime was
also performed.

This idea of combining the initiation process and the growth process is
adopted in this paper. Some of the works that use this second approach assume
that the degradation processes arrive either at a constant rate (using the
Poisson process [20]) or governed by a deterministic intensity function (using
the non-homogeneous Poisson process [5]). If the degradation processes appear
following a non-homogeneous Poisson process (NHPP) and all the degradation
processes degrade following the same degradation mechanism, the combined point
process $\left\\{S_{i}+D_{i}\right\\}$ follows a non-homogeneous Poisson
process. The non-homogeneous Poisson assumption facilitates explicit analysis
and it has been applied to real data in order to model the corroding gas
pipeline system over its lifetime [35]. However, modelling the arrival process
as a non-homogeneous Poisson process assumes that the initiation intensity is
deterministic along the system lifetime.  The motivation of our work is the
following: external shocks can accelerate the arrival of new degradation
processes to the system hence the deterministic intensity of arrivals can not
be settled. For example, in harsh environments, shocks can induce highly
dynamic loads on structures causing cracking problems [19]. In the case of
cermets, many cracks are initiated due to mechanical shocks [22]. When shocks
affect to the intensity of arrival of new degradation processes, this
phenomenon is better captured by a shot noise Cox process instead of a non-
homogeneous Poisson process. The  shot noise Cox process has been employed as
a useful tool for modeling the impact on the system lifetime of a dynamic
environment [10], [12], [13]. In this paper, we assume that the degradation
processes initiate at random times under a shot noise Cox process. From its
introduction in reliability by Lemoine and Wenocur [26] and [27], different
authors have used this process in reliability and maintenance mainly for
modeling the increment of the failure rate or an abrupt degradation increment
due to external shocks (see [17] and [31] as recent examples of the use of
shot noise process in maintenance and reliability).

Sometimes, the degradation processes present substantial variations between
them causing different degradation patterns. Different models have been
proposed to integrate the heterogeneities under the common idea that some
parameters are process-specific and different across processes [14]. An
approach to take into account the heterogeneities is to assume that the
parameters of the model follow a random variable [15]. In the case of a gamma
process, the  gamma distribution in itself is a very attractive candidate for
the distribution of the scale parameter [25]. Assuming the gamma distribution
for the scale parameter, the joint distribution of the scale parameter and the
gamma process has a closed form expression based on the Snedecor F
distribution, which allows the computation of the failure time. This gamma
random effect has been studied by different authors ([32], [40], [37] among
others). Pulcini [30] extends the gamma random effect assuming that
heterogeneity affects to the shape parameter. In this paper, unlike the
articles mentioned above, we incorporate the random effects in the gamma
process to model the process-specific heterogeneity assuming that the scale
parameter follows a uniform distribution. As far as we are concerned, there
are not many works that model the heterogeneity using a uniform distribution.

In short, in this paper, we deal with a system subject to multiple degradation
processes. We assume that the degradation processes initiate according to a
shot noise Cox process and grows according to a stationary gamma process. The
system fails when a degradation process exceeds a failure threshold. An
inspection policy is developed for this system and it is periodically
inspected.  At each inspection, the decision on whether a preventive
replacement or a corrective replacement should be taken, is performed.  The
analytical expression for the expected cost rate is obtained.

The main contributions of this paper are summarized in the following aspects.

  * 1.

A shot noise Cox process is imposed as process of initiation of degradation
processes. It extends and generalizes the combined model of initiation and
growth that assumes a non homogeneous Poisson process of arrivals. The
methodology used in this paper allows to find easy-to-evaluate expressions for
reliability expressions.

  * 2.

A random effect model is incorporated to the model with the novelty of the
probability distribution used to deal with the heterogeneity between
processes. In this paper, a uniform distribution is used to model the inverse
of the sc. It allows to evaluate the moments associated to this process and
compare the variance between models with heterogeneity and without
heterogeneity.

  * 3.

Providing the analytic cost model for the combined model of arrivals and
growth. The proposed model allows to get the analytic expressions for the
quantities related to the maintenance.

This remainder of the paper is structured as follows. In Sections 2 and 3, the
arrival and growth processes are described. In Section 4, the combined process
of initiation and growth is analyzed. Section 4 explains also the
heterogeneity model. In Section 5, the survival time is obtained and  it is
shown that the failure time distribution is increasing failure rate (IFR).
Section 6 is devoted to the maintenance analysis. Section 7 shows numerical
examples and last section conclude.

##  2 Arrival processes

We assume that a system is working in a dynamic environment and it is subject
to external shocks. The external shocks arrive to the system according to a
Poisson process with deterministic rate $\mu$. Suppose that if a shock occurs
at epoch $T_{1}$, then at time
$T_{1}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}$
the contribution of the shock to the arrival of a degradation process is
$h({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})=\exp(-\delta({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}-T_{1}))$
with $\delta>0$. Let $T_{i},\ i=1,2,\ldots$ be the arrival times of the
homogeneous Poisson process with rate parameter $\mu$ and let
$N(s)=\sum_{i=1}^{\infty}{\bf
1}_{\\{T_{i}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}\\}}$
be the counting process associated with the homogeneous Poisson process. Then,
the failure rate of the shot noise Cox process  at time $s$,
$\lambda^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})$,
is given by

$\lambda^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})=\lambda_{0}(s)+\sum_{i=1}^{N({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})}\exp{(-\delta({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}-T_{i}))},\quad
s\geq 0,$ (1)

where the deterministic function $\lambda_{0}(s)>0$ provides a Poisson base
level for the process. In absence of external shocks or if
$\delta\rightarrow\infty$, a shot noise Cox process reduces to a non
homogeneous Poisson process.

Since  the intensity $\lambda^{*}$ given by Eq. (1) is stochastic, the
expectation  at time $s$ is given by,

$\displaystyle\mathbb{E}(\lambda^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}))$
$\displaystyle=$
$\displaystyle\lambda_{0}(s)+\mathbb{E}\left(\sum_{i=1}^{N({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})}\exp{(-\delta({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}-T_{i}))}\right)$
(2) $\displaystyle=$
$\displaystyle\lambda_{0}(s)+\frac{n}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}\delta}\left(1-\exp(-\delta{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})\right),$

using that, conditional on $\left\\{N(s)=n\right\\}$,  the vector
$(T_{1},T_{2},\ldots,T_{n})$ has the same distribution as the order statistics
of sample $(U_{1},U_{2},\ldots,U_{n})$ of size $n$ from the distribution

$P(U\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t})=\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}},\quad
0\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s},$

where
$0<t_{1}<t_{2}<\ldots<t_{n}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}$
(see [27] for more details).

Next, we compute the expected number of degradation processes at time
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}$.
For that, the following lemma, valid for all counting process whose proof is
given in [33] (pp. 335-336), is used.

######  Lemma 1.

Let
$\lambda^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})$,
with
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}\geq
0$, be the random intensity function of the counting process
$\left\\{N^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}\geq
0\right\\}$ having $N^{*}(0)=0$. Then

$\displaystyle\mathbb{E}[N^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})]$
$\displaystyle=$
$\displaystyle\int_{0}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}}\mathbb{E}[\lambda^{*}(u)]~{}du.$

Using Lemma 1, the expected number of degradation processes at time
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}$
is given by

$\displaystyle\mathbb{E}[N^{*}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})]$
$\displaystyle=$
$\displaystyle\int_{0}^{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}\mathbb{E}\left[\lambda^{*}(u)\right]~{}du$
$\displaystyle=$
$\displaystyle\Lambda_{0}(s)+\frac{\mu{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s}}{\delta}+\frac{\mu}{\delta^{2}}\left(\exp{(-\delta{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}s})}-1\right).$

where

$\Lambda_{0}(s)=\int_{0}^{s}\lambda_{0}(u)du.$

##  3 The stochastic process of growth

Once a degradation process arrives to the system, the process of growth is
activated. We assume that degradation processes grow independently each other.
Due to its mathematical properties, a gamma process is used as mathematical
model of the growth. The gamma process is a stochastic process with
independent gamma-distributed increments. The gamma process with shape
function $\alpha(t)>0$ and scale parameter $\beta>0$ is a continuous-time
stochastic process $\left\\{X(t),t\geq 0\right\\}$ with the following
properties:

  1. 1.

$X(0)=0$ with probability one.

  2. 2.

$X(t_{2})-X(t_{1})\sim Gamma(\alpha(t_{2})-\alpha(t_{1}),\beta)$ for
$t_{1}\leq t_{2}$.

  3. 3.

$X(t)$ has independent increments.

Recall that a random variable $X$ has a gamma distribution with shape
parameter $\alpha>0$ and scale parameter $\beta>0$ if its probability density
function is given by

$f_{\alpha,\beta}(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}\exp^{-\beta
x},\quad x>0,$ (3)

where

$\Gamma(\alpha)=\int_{0}^{\infty}t^{\alpha-1}e^{-t}~{}dt.$

We assume in this paper that the system fails when the deterioration level of
a degradation process first exceeds a failure threshold $L$. Suppose that a
degradation process $X(t)$ starts at time 0 and it grows according to a
homogeneous gamma process with parameters $\alpha$ and $\beta$. Let
$\sigma_{L}$ be the first time at which this degradation process exceeds $L$.
The random variable $\sigma_{L}$ is known as the first hitting time
distribution and it has the following cumulative distribution function

$F_{\sigma_{L}}(t)=P(X(t)\geq L)=\int_{L}^{\infty}f_{\alpha
t,\beta}(x)~{}dx=\frac{\Gamma(\alpha t,\beta L)}{\Gamma(\alpha t)},$ (4)

for $t\geq 0$ where $f_{\alpha t,\beta}$ is given by (3) and

$\Gamma(\alpha,x)=\int_{x}^{\infty}z^{\alpha-1}e^{-z}~{}dz,$ (5)

denotes the incomplete gamma function for $x\geq 0$ and $\alpha>0$.

For subsequent analysis, the distribution of $\sigma_{L}-\sigma_{M}$ is used
in this paper for two degradation levels $M<L$. According to [7], the survival
function of this variable is given by

$\bar{F}_{\sigma_{L}-\sigma_{M}}(t)=\int_{x=0}^{\infty}\int_{y=M}^{\infty}f_{\sigma_{M},X(\sigma_{M})}(x,y)F_{\alpha
t,\beta}(L-y)~{}dy~{}dx,$ (6)

where $F_{\alpha t,\beta}$ denotes the distribution function of a gamma
distribution with parameters $\alpha t$ and $\beta$ and
$f_{\sigma_{M},X(\sigma_{M})}$ denotes the joint density function of
$(\sigma_{M},X(\sigma_{M}))$ provided in [3] as

$f_{\sigma_{M}^{h},X(\sigma_{M}^{h})}(x,y)=\int_{0}^{\infty}\mathbf{1}_{\left\\{M\leq
y<M+s\right\\}}f_{\alpha x,\beta}(y-s)\mu(ds),$

and $\mu(ds)$ denotes the Lévy measure of the gamma process with parameters
$\alpha$ and $\beta$ given by

$\mu(ds)=\alpha\frac{e^{-\beta s}}{s},\quad s>0.$

Hence,

$\displaystyle\bar{F}_{\sigma_{L}-\sigma_{M}}(t)$ $\displaystyle=$
$\displaystyle\int_{x=0}^{\infty}~{}dx\int_{s=0}^{\infty}~{}ds\int_{y=M}^{M+s}~{}dyf_{\alpha
x,\beta}(y-s)\mu(ds)F_{\alpha t,\beta}(L-y)~{}dy~{}dx$ $\displaystyle=$
$\displaystyle\int_{x=0}^{\infty}~{}dx\int_{s=0}^{\infty}~{}ds\int_{y=M}^{M+s}~{}dy\frac{\alpha(y-s)^{\alpha
x-1}}{s\Gamma(\alpha x)}\int_{u=0}^{L-y}\frac{u^{\alpha
t-1}\exp(-\beta(y+u))\beta^{\alpha x+\alpha t}}{\Gamma(\alpha t)}~{}du$

The model considered above assume that all the degradation processes are
independent and identically distributed. However, in some cases, the
degradation processes degrade at different rates even though no differences in
the environment is present [40]. In the next section, a  random effects
process-specific is used to model such variability.

###  3.1 Random effects

The random effects model is a useful tool for modeling the variability in the
different degradation rates and the heterogeneity among different degradation
processes. In this section we assume that the degradation processes can be
described by a gamma process with random effects.  To implement this random
effects model, we assume that the scale parameter of the gamma process is
random. It means that both the mean and the variance of the process are
affected by the random effect parameter.  Some authors have previously studied
the random effects model in the gamma process using different probability
distributions for the scale parameter. One of the most popular is the gamma
distribution since it provides a closed-form expression using a Fisher
distribution [25].

Inspired by the simplicity of its distribution and by the fact that in
Bayesian theory, when no a prior statistical information about the parameters
is given uniform distribution is used ([18]), a uniform distribution is used
to model the inverse of the scale parameter of the gamma process.

Let $\left\\{X_{h}(t),t\geq 0\right\\}$ be a gamma process with random effects
$\beta$ that controls the heterogeneity among the different degradation
processes. We assume that $\left\\{X_{h}(t),\,t\geq 0\right\\}$ is a
homogeneous gamma process with parameters $\alpha$ and $\beta$, where
$\beta^{-1}$ follows a uniform distribution in $(a,b)$ with $0<a<b<\infty$.
Then, the corresponding probability density function of $X_{h}(t)$ is given by

$\displaystyle f_{X_{h}(t)}(u)$ $\displaystyle=$
$\displaystyle\int_{a}^{b}f_{\alpha t,\beta}(u)\frac{1}{b-a}~{}d\beta^{-1}$
(7) $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{(b-a)}\frac{\Gamma(\alpha
t-1,u/b)-\Gamma(\alpha t-1,u/a)}{\Gamma(\alpha t)}},\quad u\geq 0,$

where $\Gamma(\cdot,\cdot)$ denotes the upper incomplete gamma function given
by Eq. (5).

Figure 1 plots the density $f_{X_{h}(t)}(u)$ for $t=5$ for a stationary gamma
process with shape parameter $\alpha=1$ and scale parameter $\beta^{-1}$
following a uniform distribution in $(1,b)$. As we can see, the kurtosis of
this distribution increases with respect to $b$.

Figure 1: Density function $f_{X_{h}(t)}$ for $\alpha=1$, $t=5$ and
$\beta^{-1}$ uniform in $(1,b)$.

Let $\sigma_{L}^{h}$ be the first time in which the process with random
effects exceeds the failure threshold $L$. The survival function of
$\sigma_{L}^{h}$ can be evaluated as

$\displaystyle{F}_{\sigma_{L}^{h}}(t)$ $\displaystyle=$ $\displaystyle
P(X_{h}(t)\geq L)$ $\displaystyle=$
$\displaystyle\int_{L}^{\infty}f_{X_{h}(t)}(u)~{}du,$

where $f_{X_{h}(t)}(u)$ is given by (7). In this random effects model, each
process degrades according to a gamma process but the degradation behavior
differs from process to process and the resulting process is not longer gamma
distributed.

Moments of $X_{h}(t)$ can be evaluated as

$\displaystyle E[X_{h}^{n}(t)]$ $\displaystyle=$
$\displaystyle\prod_{n^{*}=0}^{n-1}\left(\alpha t+n^{*}\right)E[1/\beta^{n}]$
$\displaystyle=$ $\displaystyle\prod_{n^{*}=0}^{n-1}\left(\alpha
t+n^{*}\right)\frac{1}{n+1}\sum_{k=0}^{n}a^{k}b^{n-k}.$

Hence, the expectation of $X_{h}(t)$ is given by

$\displaystyle\mathbb{E}(X_{h}(t))$ $\displaystyle=$
$\displaystyle\frac{\alpha t(a+b)}{2},$

with variance

$\displaystyle Var(X_{h}(t))$ $\displaystyle=$
$\displaystyle\mathbb{E}(X_{h}(t)^{2})-(\mathbb{E}(X_{h}(t)))^{2}$
$\displaystyle=$ $\displaystyle\frac{(\alpha t+(\alpha
t)^{2})(b^{2}+ab+a^{2})}{3}-\frac{(\alpha t)^{2}(a+b)^{2}}{4}$
$\displaystyle=$ $\displaystyle\frac{\alpha
t(b^{2}+ab+a^{2})}{3}+\frac{(\alpha t)^{2}}{12}(a-b)^{2}.$

For fixed $t<\infty$, and assuming that $0<a<b<\infty$, $X_{h}(t)$ has finite
mean and variance.

Expectation and variance of the random effects process
$\left\\{X_{h}(t),\,t\geq 0\right\\}$ increase with the time. As in the case
without heterogeneities, expectation increases linearly with the time. In the
case of the variance, unlike the case without heterogeneities, variance does
not increase linearly with the time. Figures 2a and 2b show the expectation
and variance of the process $\left\\{X_{h}(t),\,t\geq 0\right\\}$ for a
stationary gamma process with shape parameter $\alpha=1$ and $\beta^{-1}$
following a uniform distribution in $(1,b)$.

(a) Expectation

(b) Variance

Figure 2: Expectation and variance of $X_{h}(t)$ versus $t$ for $\alpha=1$ and
$\beta^{-1}$ uniform in $(1,b)$

It is well known that, in a gamma process $X(t)$ without heterogeneity, the
ratio $\text{Var}(X(t))/\mathbb{E}(X(t))$ does not depend on time. However, in
the gamma process with heterogeneity presented in this paper, the quotient

$\displaystyle\frac{\text{Var}(X_{h}(t))}{\mathbb{E}(X_{h}(t))}$
$\displaystyle=$ $\displaystyle\frac{2(1+\alpha
t)(b^{2}+ab+a^{2})}{3(a+b)}-\frac{\alpha t(a+b)}{2}$ (8) $\displaystyle=$
$\displaystyle\frac{4(b^{2}+ab+a^{2})+\alpha t(b-a)^{2}}{6(a+b)},$

is no longer constant (increases with $t$). As [30] claims, the presence of a
noticeable heterogeneity should ensure that the ratio (8) increases with the
time. Figure 3 shows the ratio given by Eq. (8) of the process
$\left\\{X_{h}(t),\,t\geq 0\right\\}$ for a stationary gamma process with
shape parameter $\alpha=1$ and $\beta^{-1}$ following a uniform distribution
in $(1,b)$. As we can see, the ratio is increasing with the time.

Figure 3: $Var(X_{h}(t))/E(X_{h}(t))$ versus $t$ with $\beta^{-1}$ uniform in
$(1,b)$

The random effects model with heterogeneity modelled using a uniform
distribution is next compared to a random effects model with heterogeneity
modelled using a gamma distribution and compared to a model without
heterogeneity. For that, let $\left\\{X_{h}^{(1)}(t),\,t\geq 0\right\\}$,
$\left\\{X_{h}^{(2)}(t),\,t\geq 0\right\\}$ and $\left\\{X(t),\,t\geq
0\right\\}$ be three stochastic processes where $X_{h}^{(1)}(t)$ denotes the
degradation level at time $t$ for a gamma process with shape parameter
$\alpha$ and scale parameter $\beta_{1}$ where $\beta_{1}^{-1}$ follows a
uniform distribution $U(a,b)$, $X_{h}^{(2)}(t)$ denotes the degradation level
at time $t$ for a gamma process with shape parameter $\alpha$ and scale
parameter $\beta_{2}$ where $\beta_{2}$ follows a gamma distribution with
parameter $k_{1}$ and $k_{2}$ $Gamma(k_{1},k_{2})$ and $X(t)$ denotes the
degradation level at time $t$ for a gamma process with shape parameter
$\alpha$ and scale parameter $\beta_{3}$ with deterministic $\beta_{3}$. To
compare the three processes, firstly the expectations are standardized
imposing

$E[X_{h}^{(1)}(t)]=E[X_{h}^{(2)}(t)]=E[X(t)],\quad t\geq 0.$

To get the same expectation, the parameters of the three processes fulfill

$\frac{1}{\beta_{3}}=\frac{a+b}{2}=\frac{k_{1}}{k_{2}}.$ (9)

Since the three processes have the same expectation, we get that the
comparison in terms of variances is reduced to the comparison of the moments
of second order $\mathbb{E}[X_{h}^{(1)}(t)^{2}]$,
$\mathbb{E}[X_{h}^{(2)}(t)^{2}]$ and $\mathbb{E}[X(t)^{2}]$. We get that

$\displaystyle\mathbb{E}[X_{h}^{(1)}(t)^{2}]$ $\displaystyle=$
$\displaystyle\left(\alpha t+(\alpha
t)^{2}\right)\mathbb{E}\left[\frac{1}{\beta_{1}^{2}}\right]=\left(\alpha
t+(\alpha t)^{2}\right)\frac{b^{2}+ab+a^{2}}{3}$
$\displaystyle\mathbb{E}[X_{h}^{(2)}(t)^{2}]$ $\displaystyle=$
$\displaystyle\left(\alpha t+(\alpha
t)^{2}\right)\mathbb{E}\left[\frac{1}{\beta_{2}^{2}}\right]=\left(\alpha
t+(\alpha t)^{2}\right)\frac{k_{1}(k_{1}+1)}{k_{2}}$
$\displaystyle\mathbb{E}[X(t)^{2}]$ $\displaystyle=$
$\displaystyle\left(\alpha t+(\alpha
t)^{2}\right)\mathbb{E}\left[\frac{1}{\beta_{3}^{2}}\right]$

Comparing $X_{h}^{(1)}(t)$ and $X_{h}^{(2)}(t)$ with parameters fulfilling Eq.
(9) in terms of variances, we get that

$Var(X_{h}^{(1)}(t))\geq(\leq)\,Var(X_{h}^{(2)}(t))\Leftrightarrow\frac{2(b^{2}+ab+a^{2})}{3(a+b)}-1\geq(\leq)\,k_{1}.$

Hence, the variances of the two processes $Var(X_{h}^{(1)}(t))$ and
$Var(X_{h}^{(2)}(t))$ depend on the parameters of the processes. Comparing the
random effect model with uniform distribution and the process without
heterogeneity, we get that

$\displaystyle\mathbb{E}[X_{h}^{(1)}(t)^{2}]$ $\displaystyle=$
$\displaystyle\left(\alpha t+(\alpha t)^{2}\right)\frac{b^{2}+ab+a^{2}}{3}$
$\displaystyle\mathbb{E}[X(t)^{2}]$ $\displaystyle=$
$\displaystyle\left(\alpha t+(\alpha t)^{2}\right)\frac{(a+b)^{2}}{4}$

hence $Var(X_{h}^{(1}(t))\geq Var(X(t))$ for all $t$ and parameters fulfilling
Eq. (9). Hence for two processes with the same expectation, a gamma process
without heterogeneity shows a lower variance with respect to a gamma process
with random effects model under a uniform distribution.

Expression (7) allows us to obtain the likelihood function for observed data.
Suppose now that $X_{h}(t)$ is recorded at times $t_{1},t_{2},\ldots,t_{n}$
with observations $x_{1},x_{2},\ldots,x_{n}$. The $x$-increments are defined
by $\Delta x_{j}=x_{j}-x_{j-1}$ for $j=1,2,\ldots,n$ and $x_{0}=0$ and the
$t$-increments are defined by $\Delta t_{j}=t_{j}-t_{j-1}$ with $t_{0}=0$.
Using (7) for fixed $\alpha$ and for a realization of $\beta^{-1}$, the joint
density of $\left(\Delta x_{1},\Delta x_{2},\ldots,\Delta x_{n}\right)$ is
given by the product

$f(\Delta x_{1},\Delta x_{2},\ldots,\Delta x_{n})=\frac{\prod_{j=1}^{n}(\Delta
x_{j})^{\alpha\Delta t_{j}-1}}{\prod_{j=1}^{n}\Gamma(\alpha\Delta
t_{j})}\beta^{\alpha t_{n}}\exp(-\beta x_{n}).$ (10)

Suppose that $m$ degradation processes $\left\\{X_{h,i}(t),t\geq 0\right\\}$
for $i=1,2,\ldots,m$ are observed with the $i$-th process observed at times
$t_{0,i}<t_{1,i},\ldots,<t_{n_{i},i}$ with values
$x_{0,i}<x_{1,i},\ldots,<x_{n_{i},i}$. Analogously to [25], the joint density
of the increments $\Delta x_{j,i}=x_{j,i}-x_{j-1,i}$ for the $i$ process is
given by (10) replacing $\Delta x_{j}$ by $\Delta x_{j,i}$, $\Delta t_{j}$ by
$\Delta t_{j,i}$. The product across $i=1,2,\ldots,m$ and the integration with
respect to $\beta$ gives the likelihood function of the observed data in the
$m$ processes.

$L(\alpha,a,b)=\frac{1}{(b-a)^{m}}\prod_{i=1}^{m}\frac{\prod_{j=1}^{n_{i}}(\Delta
x_{j,i})^{\alpha\Delta t_{j,i}-1}}{\prod_{j=1}^{n_{i}}\Gamma(\alpha\Delta
t_{j,i})}\frac{\Gamma(\alpha t_{n,i}-1,x_{n,i}/b)-\Gamma(\alpha
t_{n,i}-1,x_{n,i}/a)}{x_{n,i}^{\alpha t_{n,i}}}.$ (11)

Inference about $(\alpha,a,b)$ can be based on the usual maximum likelihood
procedures using Eq.(11).

######  Example 1.

Different degradation processes have been simulated from $t=0$ to $t=30$ for
$\alpha=1.5$ and $\beta^{-1}$ uniform in $(1-0.3,1+0.3)$. The trajectories are
plotted in Figures 4a and 5a for six degradation processes and twenty six
degradation processes respectively.

For fixed $\alpha=1.5$ and $\beta^{-1}$ uniform in
$(1-\alpha^{*},1+\alpha^{*})$, Figures 4b and 5b show the negative of the
logarithm of the likelihood of the six simulation processes versus
$\alpha^{*}$. As we can check, the minimum of this function is reached for
$\alpha^{*}=0.39$ with a value of -473.5824 (for 6 simulated processes) and
$\alpha^{*}=0.302$ with a value of the negative of the likelihood of -431.5405
for 26 degradation processes.

(a) Trajectories

(b) Negative of the log likelihood

Figure 4: Trajectories and log likelihood

(a) Trajectories

(b) Negative of the log likelihood

Figure 5: Trajectories and log likelihood

##  4 Combined process of arrival and growth

Since the shot noise process with intensity given by Eq. (1) is a Cox process,
this section shows some previous well-known results about general Cox
processes. These results are used to evaluate quantities of interest of the
model.  The degradation processes start at random times $S_{1},S_{2},\ldots$
following a shot noise Cox process with intensity given by (1) and then they
grow independently following a homogeneous gamma process with parameters
$\alpha$ and $\beta$. Next, we shall analyse the combined process of arrival
and growth.

We define the deterioration level of the degradation process $k$ at time $t$
as

$X_{k}(t)=X^{(k)}(t-S_{k}),\quad t\geq S_{k},$

where $X^{(k)}$ denotes a homogeneous gamma process with parameters $\alpha$
and $\beta$. We set $X^{(i)}$ with $i\in N$ to be i.i.d of a stationary gamma
process with parameters $\alpha$ and $\beta$.

Let $W_{k}$ be the time where the deterioration level of the degradation
process $k$ exceeds the threshold $L$. Hence, $W_{k}$ is given by

$W_{k}=S_{k}+\sigma_{L},\ \ \text{for}\ k=1,2,\ldots$

Figure 6 shows a realization of the combined process.

Figure 6: Realization of the combined process.

Let $N_{L}(t)$ be the number of degradation processes whose deterioration
levels exceed the failure threshold $L$ at time $t$, that is,

$N_{L}(t)=\sum_{k=1}^{\infty}\mathbf{1}_{\left\\{W_{k}\leq t\right\\}}.$ (12)

The distribution of the counting process $\left\\{N_{L}(t),\,t\geq 0\right\\}$
is next obtained.

Before obtaining the distribution of $\left\\{N_{L}(t),\,t\geq 0\right\\}$,
some definitions and properties of the Cox processes are first recalled. These
properties allow to obtain the combined expression for arrival and growth of
the degradation processes.

Because Cox process are essentially Poisson processes, most results for
Poisson processes have counterparts for Cox processes using an expectation to
take into account the stochastic intensity. The preservation of basic point
process operations is an example of these results. We focus on the random
translation of a point process, the so-called displaced process.

######  Definition 1.

Let $\left\\{N^{*}(t),t\geq 0\right\\}$ be a counting process with occurrence
times $S_{1},S_{2},\ldots$ and suppose that
$\left\\{D_{1},D_{2},D_{3},\ldots\right\\},\quad i=1,2,\ldots$ is a sequence
of independent and identically distributed, non-negative random variables
independent of $S_{i}$. The point process $\left\\{N(t),t\geq 0\right\\}$
whose occurrence times are $\left\\{S_{i}+D_{i}\right\\}$ for $i=1,2,\ldots$
is called a displaced process.

It is well known that if $\left\\{N^{*}(t),t\geq 0\right\\}$ follows a non-
homogeneous Poisson process with intensity $\lambda^{*}(t)$, then the
displaced process $\left\\{N(t),t\geq 0\right\\}$ with occurrence times given
by $\left\\{T_{i}+D_{i}\right\\}$ is also a non homogeneous Poisson process
[5] with intensity

$\lambda(t)=\int_{0}^{t}\lambda^{*}(u)f(t-u)~{}du,\quad t\geq 0,$

where $f$ denotes the common density of the random variables $D_{i}$.

This preservation of the random translation of a non homogeneus Poisson
process has the counterpart for the Cox process. It is explained in the
following result (see References [34]and [29] for more details).

######  Lemma 2.

If $\left\\{N^{*}t),t\geq 0\right\\}$ is a Cox process with stochastic
intensity $\lambda^{*}(t)$, then the displaced process $\left\\{N(t),t\geq
0\right\\}$ with occurrence times given by $\left\\{T_{i}+D_{i}\right\\}$ is
also a Cox process with stochastic intensity

$\lambda(t)=\int_{0}^{t}\lambda^{*}(u)f(t-u)~{}du,\quad t\geq 0,$

where $f$ denotes the common density of the random variables $D_{i}$.

Since, in this paper, the arrival process (the shot noise Cox process) is a
Cox process, using Lemma 2, the mathematical model for the counting process
$\left\\{N_{L}(t),\,t\geq 0\right\\}$ given by Eq. (12) is obtained and it is
shown in the next result.

######  Lemma 3.

Let $S_{1},S_{2},\ldots$ be the occurrence times of the shot noise Cox process
with stochastic intensity given by (1), then the displaced process
$\left\\{N_{L}(t),t\geq 0\right\\}$ with occurrence times given by
$\left\\{S_{i}+\sigma_{L}\right\\}$ is also a Cox process with stochastic
intensity

$\lambda_{L}(t)=\int_{0}^{t}\lambda^{*}(u)f_{\sigma_{L}}(t-u)~{}du,\quad t\geq
0,$ (13)

where $f_{\sigma_{L}}$ denotes the density of the first hitting time
$\sigma_{L}$ with distribution given by Eq. (4).

The expectation of the stochastic intensity $\lambda_{L}(t)$ of the combined
process of arrivals and growth is obtained.

######  Lemma 4.

The expectation of the stochastic intensity of the displaced process
$\left\\{N_{L}(t),t\geq 0\right\\}$ with occurrence times given by
$\left\\{S_{i}+\sigma_{L}\right\\}$ and where $S_{i}$ are given by a shot
noise Cox process with intensity given by Eq. (1) is given by

$\mathbb{E}[\lambda_{L}(t)]=\lambda_{0}F_{\sigma_{L}}(t)+\mu\int_{0}^{t}H(u)f_{\sigma_{L}}(t-u)~{}du,$

where $f_{\sigma_{L}}$ is the density function of the hitting time
$\sigma_{L}$ and $H(t)$ is given by

$H(t)=\int_{0}^{t}e^{-\delta u}~{}du.$

###### Proof.

As we explained before, $\left\\{N_{L}(t),t\geq 0\right\\}$ is a Cox process
with intensity $\lambda_{L}(t)$ given by (13) with expectation

$\mathbb{E}[\lambda_{L}(t)]=\int_{0}^{t}\mathbb{E}[\lambda^{*}(u)]f_{\sigma_{L}}(t-u)~{}du,$

and, using $\mathbb{E}[\lambda^{*}(u)]$ given in (2), the result is obtained
using Fubini theorem since $\lambda^{*}$ is integrable. ∎

Finally, Lemma 4 allows us to obtain the expected number of degradation
processes that exceed the failure threshold $L$ at time $t$ given by Eq. (12).
To derive $\mathbb{E}[N_{L}(t)]$, Lemma 1 is used.

$\displaystyle\mathbf{E}\left[N_{L}(t)\right]=\int_{0}^{t}\mathbb{E}[\lambda_{L}(s)]~{}ds=\int_{0}^{t}\mathbb{E}[\lambda^{*}(u)]F_{\sigma_{L}}(t-u)~{}du,$

where $\mathbb{E}[\lambda^{*}(u)]$ is given by (2).

Notice that Lemmas 3 and 4 have a counterpart in the random effects case
replacing $f_{\sigma_{L}}$ by $f_{\sigma_{L}^{h}}$ which is the density
probability function of $\sigma_{L}^{h}$.

Next section analyzes the failure time for this system. The failure time
distribution is obtained and it is shown that it is increasing failure rate.

##  5 Time to the system failure

The system can be considered to be failed if at least one of the degradation
processes reaches its failure threshold. We assume that the failure threshold
is $L$ for all the degradation processes.

Let $W_{[1]}$ be the instant at which, for the first time, the deterioration
level of a degradation process exceeds the  failure threshold L,

$W_{[1]}=\min_{i=1,2,\ldots}\left\\{W_{i}\right\\}.$

The probability distribution of $W_{[1]}$ is next obtained.

######  Lemma 5.

Let $W_{[1]}$ be the instant at which, for the first time, the deterioration
level of a degradation process exceeds the failure threshold $L$. Then, the
survival function of $W_{[1]}$ is given by

$\bar{F}_{W_{[1]}}(t)=\exp\left(-\lambda_{0}\int_{0}^{t}F_{\sigma_{L}}(u)~{}du-\mu\int_{0}^{t}\left(1-\exp\left(-\int_{0}^{u}e^{-\delta
w}F_{\sigma_{L}}(u-w)~{}dw\right)\right)~{}du\right)$

where $F_{\sigma_{L}}$ is given by (4).

###### Proof.

As we obtained above, counting process $\left\\{N_{L}(t),t\geq 0\right\\}$ is
a Cox process with intensity $\lambda_{L}(t)$ given by (13). Hence

$\displaystyle\bar{F}_{W_{[1]}}(t)$ $\displaystyle=$ $\displaystyle
P(N_{L}(t)=0)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\exp\left\\{-\int_{0}^{t}\lambda_{L}(u)~{}du\right\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\exp\left\\{-\int_{0}^{t}\lambda^{*}(u)F_{\sigma_{L}}(t-u)~{}du\right\\}\right]$
$\displaystyle=$
$\displaystyle\exp\left\\{-\lambda_{0}\int_{0}^{t}F_{\sigma_{L}}(u)~{}du\right\\}\mathbb{E}\left[\exp\left\\{-\int_{0}^{t}\sum_{i=1}^{N(t)}e^{-\delta(u-T_{i})}F_{\sigma_{L}}(t-u)\mathbf{1}_{\left\\{u>T_{i}\right\\}}~{}du\right\\}\right].$

Given $N(t)=n$ with $n>0$, we get that

$\displaystyle P\left(W_{[1]}>t|N(t)=n\right)$ $\displaystyle=$ $\displaystyle
C_{1}(t)\mathbb{E}\left[\prod_{i=1}^{n}\exp\left\\{-\int_{0}^{t}e^{-\delta(u-T_{i})}F_{\sigma_{L}}(t-u)\mathbf{1}_{\left\\{u>T_{i}\right\\}}~{}du\right\\}\right]$
$\displaystyle=$ $\displaystyle
C_{1}(t)\mathbb{E}\left[\left(\exp\left\\{-\int_{0}^{t}e^{-\delta(u-U)}F_{\sigma_{L}}(t-u)\mathbf{1}_{\left\\{u>U\right\\}}~{}du\right\\}\right)\right]^{n},$

where

$C_{1}(t)=\exp\left(-\lambda_{0}\int_{0}^{t}F_{\sigma_{L}}(u)~{}du\right),\quad
t\geq 0,$ (14)

and $U$ is a uniform variable in the interval $(0,t)$. Then

$\displaystyle P\left(W_{[1]}>t|N(t)=n\right)$ $\displaystyle=$ $\displaystyle
C_{1}(t)\left(\int_{s=0}^{t}\frac{1}{t}\exp\left\\{-\int_{s}^{t}e^{-\delta(u-s)}F_{\sigma_{L}}(t-u)~{}du\right\\}~{}ds\right)^{n},$

therefore

$\displaystyle\bar{F}_{W_{[1]}}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}P(W_{[1]}>t|N(t)=n)P(N(t=n))$
$\displaystyle=$ $\displaystyle C_{1}(t)\exp(-\mu
t)+\sum_{n=1}^{\infty}P(W_{[1]}>t)P(N(t=n))$ $\displaystyle=$ $\displaystyle
C_{1}(t)\exp(-\mu t)+C_{1}(t)\exp(-\mu
t)\sum_{n=1}^{\infty}\frac{A(t)^{n}}{t^{n}}\frac{(\mu t)^{n}}{n!}$

where $A(t)$ is given by

$A(t)=\int_{0}^{t}\exp\left(-\int_{0}^{x}e^{-\delta
w}F_{\sigma_{L}}(x-w)~{}dw\right)~{}dx.$

Hence

$\displaystyle\bar{F}_{W_{[1]}}(t)$ $\displaystyle=$ $\displaystyle
C_{1}(t)\exp(-\mu t)+C_{1}(t)\exp(-\mu t)\left(\exp(A(t)\mu)-1\right)$
$\displaystyle=$ $\displaystyle C_{1}(t)\exp(-\mu t))\exp(A(t)\mu$
$\displaystyle=$ $\displaystyle C_{1}(t)C_{2}(t)$

where $C_{1}(t)$ is given by (14) and $C_{2}(t)$ equals to

$C_{2}(t)=\exp\left(-\mu\int_{0}^{t}\left(1-\exp\left(-\int_{0}^{x}e^{-\delta
w}F_{\sigma_{L}}(x-w)~{}dw\right)\right)~{}dx\right),$ (15)

and the result holds. ∎

######  Lemma 6.

Let $W_{[1]}$ be the instant at which, for the first time, the deterioration
level of a degradation process exceeds the failure threshold $L$. Then
$W_{[1]}$ is increasing failure rate (IFR).

###### Proof.

Using Lemma 5, the failure rate function of  $W_{[1]}$ is given by

$r_{W_{[1]}}(t)=\frac{-C_{1}^{\prime}(t)}{C_{1}(t)}+\frac{-C_{2}^{\prime}(t)}{C_{2}(t)},$

where $C_{1}$ and $C_{2}$ are given by (14) and (15) respectively.  By making
certain calculations, we get that

$r_{W_{[1]}}(t)=\lambda_{0}F_{\sigma_{L}}(t)+\mu\left(1-\exp\left(-\int_{0}^{t}e^{-\delta
w}F_{\sigma_{L}}(t-w)~{}dw\right)\right),$

with derivative

$r_{W_{[1]}}^{\prime}(t)=\lambda_{0}f_{\sigma_{L}}(t)+\mu\exp\left(-\int_{0}^{t}\exp^{-\delta
w}F_{\sigma_{L}}(t-w)~{}dw\right)\left(\int_{0}^{t}e^{-\delta
w}f_{\sigma_{L}}(t-w)~{}dw\right)$

and it is non-decreasing in $t$ independent of the monotony of the failure
rate of $\sigma_{L}$. ∎

Notice that Lemmas 5 and 6 have a counterpart replacing $F_{\sigma_{L}}$ by
$F_{\sigma_{L}^{h}}$.

The limit of $r_{W_{[1]}}(t)$ is given by

$\lim_{t\rightarrow\infty}r_{W_{[1]}}(t)=\lambda_{0}+\mu\left(1-\exp(-1/\delta)\right).$

The increasing failure rate implies that the preventive maintenance is
potentially worth implementing to improve the system reliability.

##  6 Maintenance policy

It is assumed that the degradation of the degradation processes cannot be
directly and instantaneously observed. The system is monitored through
periodic inspections which reveal the exact state of the degradation processes
without error.  Failures are non self-announcing, so that system failures are
only detected through inspections. This paper makes the following assumptions.

  1. 1.

The system is inspected each $T$ time units and the following decisions are
taken.

     * (a)

If the degradation level of a degradation process exceeds a preventive
maintenance threshold $M$ but the system is not failed, a preventive
replacement is performed and the system is replaced by a new one.

     * (b)

If the system is failed, a corrective replacement is performed and the system
is replaced by a new one.

     * (c)

If there is no degradation processes in the system or their degradation levels
do not exceed the preventive maintenance threshold, the system is left as it
is.

  2. 2.

Replacements are immediately performed and the duration time of these
replacements is negligible.

  3. 3.

A sequence of costs is associated with the different maintenance actions:

     * (a)

A cost of $C_{c}$ monetary units is incurred when a corrective replacement is
performed.

     * (b)

A cost of $C_{p}$ monetary units is incurred when a preventive replacement is
performed.

     * (c)

The cost of the inspections is equal to $C_{I}$ monetary units.

     * (d)

The downtime cost is equal to $C_{d}$ monetary units per  time units.

First, let $\left\\{N_{M}(t),t\geq 0\right\\}$ be the number of degradation
processes that exceed the preventive threshold $M$ at time $t$. That is,

$N_{M}(t)=\sum_{i=1}^{\infty}\mathbf{1}_{\left\\{S_{i}+\sigma_{M}\leq
t\right\\}}.$

Using a similar reasoning to the development of $\left\\{N_{L}(t),t\geq
0\right\\}$, the point process $\left\\{N_{M}(t),t\geq 0\right\\}$ follows a
Cox process with stochastic intensity

$\lambda_{M}(t)=\int_{0}^{t}\lambda^{*}(u)f_{\sigma_{M}}(t-u)~{}du,\quad t\geq
0,$ (16)

where $\lambda^{*}(u)$ is given by (1)

Let $V_{[1]}$ be the instant at which, for the first time, the deterioration
level of a degradation process exceeds the preventive threshold $M$. Using a
similar reasoning to the method developed above, the distribution of $V_{[1]}$
can be expressed as

$\displaystyle\bar{F}_{V_{[1]}}(t)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\exp\left\\{-\int_{0}^{t}F_{\sigma_{M}}(t-u)\lambda^{*}(u)~{}du\right\\}\right],\quad
t\geq 0,$

where $F_{\sigma_{M}}$ is given by Eq. (4) replacing $L$ by $M$.

Let $R$ be the time to a replacement, that is, the time to a preventive
replacement or a corrective maintenance. Hence,

$R=\sum_{k=0}^{\infty}(k+1)T\mathbf{1}_{\left\\{kT\leq
V_{[1]}\leq(k+1)T\right\\}},$

with expectation

$\displaystyle\mathbb{E}[R]$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}(k+1)TP(kT\leq V_{[1]}\leq(k+1)T)$
$\displaystyle=$ $\displaystyle T\sum_{i=0}^{\infty}\bar{F}_{V_{[1]}}(iT).$

Let $N_{I}(R)$ be the total number of inspections performed on the system in a
replacement cycle, then

$\mathbb{E}[N_{I}(R)]=\frac{\mathbb{E}[R]}{T}=\sum_{i=0}^{\infty}\bar{F}_{V_{[1]}}(iT).$
(18)

### Preventive maintenance probability

A preventive replacement is performed at time $(k+1)T$ if the system is not
failed at time $(k+1)T$ and the preventive threshold $M$ is exceeded for the
first time in $(kT,(k+1)T]$, that is,

$\left\\{kT\leq V_{[1]}\leq(k+1)T<W_{[1]}\right\\}.$

Let $\left\\{N_{L-M}(t),t\geq 0\right\\}$ be the following counting process

$N_{L-M}(t)=\sum_{n=0}^{\infty}\mathbf{1}_{\left\\{V_{i}+(\sigma_{L}-\sigma_{M})\leq
t\right\\}},$ (19)

where the distribution of $\sigma_{L}-\sigma_{M}$ is given by Eq. (6). Since
$\left\\{V_{i}\right\\}$, for $i=1,2,\ldots$ is a Cox process with intensity
$\lambda_{M}$ given by (16), the translated process $N_{L-M}(t)$ given in (19)
is a Cox process with intensity

$\lambda_{L-M}(t)=\int_{0}^{t}\lambda_{M}(u)f_{\sigma_{L}-\sigma_{M}}(t-u)~{}du,$

where $f_{\sigma_{L}-\sigma_{M}}(\cdot)$ denotes the density of (6).

Let $P_{p}((k+1)T)$ be the preventive maintenance probability at time
$(k+1)T$. That is,

$P_{p}((k+1)T)=P\left(kT\leq V_{[1]}<(k+1)T\leq W_{[1]}\right).$

This probability is developed as follows.

$\displaystyle P_{p}((k+1)T)$ $\displaystyle=$
$\displaystyle\int_{kT}^{(k+1)T}f_{V_{[1]}}(u)P\left((k+1)T\leq
W_{[1]}|V_{[1]}=u\right)~{}du$ $\displaystyle=$
$\displaystyle\int_{kT}^{(k+1)T}f_{V_{[1]}}(u)\bar{F}_{\sigma_{L}-\sigma_{M}}((k+1)T-u)P\left(N_{L-M}(u,(k+1)T)=0\right)~{}du.$

Since $\left\\{N_{L-M}(t),\,t\geq 0\right\\}$ is a Cox process, we get that

$P(N_{L-M}(u,(k+1)T)=0)=\mathbb{E}\left[\exp\left\\{-\int_{u}^{(k+1)T}\lambda_{M}(v)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)~{}dv\right\\}\right].$

Let $g(u,(K+1)T)$ be the following expectation

$g(u,(K+1)T)=\mathbb{E}\left[\exp\left\\{-\int_{u}^{(k+1)T}\lambda_{M}(v)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)~{}dv\right\\}\right].$

Using Eq. (16)

$\lambda_{M}(t)=\lambda_{0}F_{\sigma_{M}}(t)+\sum_{i=1}^{N(t)}\int_{0}^{t}e^{-\delta(u-T_{i})}f_{\sigma_{M}}(t-u)~{}du,$

and given $N((k+1)T)=n$, for $n>0$, we get that

$\displaystyle
g_{n}(u,(K+1)T)=\exp\left\\{-\lambda_{0}\int_{u}^{(k+1)T}F_{\sigma_{M}}(v)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)~{}dv\right\\}$
$\displaystyle\mathbb{E}\left[\exp\left\\{-\int_{u}^{(k+1)T}\left(\int_{w=0}^{v}e^{-\delta(w-U)}f_{\sigma_{M}}(v-w)\mathbf{1}_{\left\\{w>U\right\\}}\right)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)\right\\}\right]^{n}$
$\displaystyle=D_{1}(u,(K+1)T)D_{2}(u,(K+1)T),$

where $U$ is a uniform variable in $(0,(k+1)T)$ and $D_{1}((K+1)T)$ and
$D_{2}((K+1)T)$ are given by

$D_{1}(u,(k+1)T)=\exp\left\\{-\lambda_{0}\int_{u}^{(k+1)T}F_{\sigma_{M}}(v)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)~{}dv\right\\},$
(20)

and

$\displaystyle D_{2}(u,(K+1)T)$ $\displaystyle=$
$\displaystyle\left(\int_{s=0}^{(k+1)T}\frac{1}{(k+1)T}\exp\left\\{-\int_{v=u}^{(k+1)T}\mathcal{L}\sigma_{M}(v-s)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)\mathbf{1}_{\left\\{v>s\right\\}}~{}dv\right\\}\right)^{n}$
$\displaystyle=$ $\displaystyle\frac{B(u,(K+1)T)^{n}}{(k+1)^{n}T^{n}}$

where

$B(u,(K+1)T)=\int_{s=0}^{(k+1)T}\exp\left\\{-\int_{u}^{(k+1)T}\mathcal{L}\sigma_{M}(v-s)F_{\sigma_{L}-\sigma_{M}}((k+1)T-v)~{}dv\right\\}~{}ds,$

and

$\mathcal{L}\sigma_{M}(x)=\int_{0}^{x}e^{-\delta z}f_{\sigma_{M}}(x-z)~{}dz.$

Finally, the preventive maintenance probability at time $(k+1)T$ is given by

$\displaystyle g(u,(K+1)T)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}g_{n}(u,(k+1)T)P[N((k+1)T)=n]$ (21)
$\displaystyle=$ $\displaystyle\exp(-\mu(K+1)T)D_{1}(u,(k+1)T)$
$\displaystyle+$
$\displaystyle\exp(-\mu(K+1)T)D_{1}(u,(k+1)T)\sum_{n=1}^{\infty}\frac{B(u,(K+1)T)^{n}}{k+1)^{n}T^{n}}\frac{\mu^{n}(k+1)^{n}T^{n}}{n!}$
$\displaystyle=$ $\displaystyle\exp(-\mu(K+1)T)D_{1}(u,(k+1)T)\exp(\mu
B(u,(k+1)T)),$

where $D_{1}(u,(k+1)T)$ is given by (20). Finally, the probability of a
maintenance preventive at time $(k+1)T$ is given by

$P_{p}((k+1)T)=\int_{kT}^{(k+1)T}f_{V_{[1]}}(u)\bar{F}_{\sigma_{L}-\sigma_{M}}((k+1)T-u)g(u,(K+1)T)~{}du.$
(22)

### Corrective maintenance probability

Let $P_{c}((k+1)T)$ be the probability of a corrective maintenance. This
probability is given by

$P_{c}((K+1)T)=P(kT\leq V_{[1]}\leq W_{[1]}\leq(K+1)T).$

Using (22), we get that

$\displaystyle
P_{c}((k+1)T)=\int_{kT}^{(k+1)T}f_{V_{[1]}}(u)\left(\int_{u}^{(K+1)T}\frac{d}{dv}\left(\bar{F}_{\sigma_{L}-\sigma_{M}}((k+1)T-v)g(v,(K+1)T)\right)~{}dv\right)~{}du,$
(23)

where $g(v,(K+1)T)$ is given by (21).

###  6.1 Expected downtime in a replacement cycle

Let $E_{d}(kT,(k+1)T)$ be the expected downtime in $(KT,(k+1)T)$. This
expectation is given by

$\displaystyle E_{d}(kT,(k+1)T)=$ (24)
$\displaystyle\int_{kT}^{(k+1)T}f_{V_{[1]}}(u)\left(\int_{u}^{(K+1)T}\frac{d}{dv}\left(\bar{F}_{\sigma_{L}-\sigma_{M}}((k+1)T-v)g(v,(K+1)T)\right)\left((K+1)T-v\right)~{}dv\right)~{}du,$

where $g(v,(K+1)T)$ is given by (21)

To evaluate the performance of the maintenance strategy, the minimization of
the asymptotic cost rate is applied. The renewal theorem (see Tijms [36] for
more details) provides a simple expression of this asymptotic cost rate since
it reduces the renewal process to the first renewal. In this paper, a renewal
cycle is defined as the period of time between two consecutive replacements of
the system. The expected cost rate for this maintenance model is given by

$C(T,M)=\frac{\mathbb{E}[C]}{\mathbb{E}[R]}$

where $\mathbb{E}[C]$ denotes the expected cost in a replacement cycle and
$\mathbb{E}[R]$ denotes the expected time to a replacement cycle. Hence

$C(T,M)=\frac{C_{c}\sum_{k=1}^{\infty}P_{c}(kT)+C_{p}\sum_{k=1}^{\infty}P_{p}(kT)+C_{I}\mathbb{E}[N_{I}]+\sum_{k=0}^{\infty}C_{d}\mathbb{E}_{d}(kT,(k+1)T)}{\mathbb{E}[R]},$
(25)

where $\mathbb{E}[R]$, $\mathbb{E}[N_{I}]$, $P_{p}(kT)$, $P_{c}(kT)$, and
$E_{d}(kT,(k+1)T)$ are given by Eq. (6), Eq. (18), Eq. (22), Eq. (23) and Eq.
(24) respectively. The search of the optimal maintenance strategy is reduced
to find the values $T_{opt}$ and $M_{opt}$ that minimize the function $C(T,M)$
given by Eq. (25). Next section focuses on numerical examples.

##  7 Numerical examples

In this section, numerical examples obtained by simulation of the stochastic
processes involved in the model are shown. The expression (25) for the
objective cost function is difficult to manage it since it evolves infinity
sums. A classical way to evaluate the objective cost function is to perform
simulations on the points of the mesh and to find the optimal combination by
visualization (see [21]) or using optimization metaheuristic methods such as
genetic algorithms ([2]). We analyze the optimal maintenance strategy
considering a model without heterogeneity and a model with heterogeneity.

###  7.1 $\beta$ deterministic

We assume that degradation processes start according to a stochastic intensity
given by

$\lambda(t)=1+\sum_{i=1}^{N(t)}e^{-0.5(t-T_{i})},\quad t\geq 0,$ (26)

where $N(t)$ denotes the Poisson process  of shocks with parameter $\mu=2$
shocks per time unit. The degradation processes degrade according to a gamma
process with shape parameter $\alpha=1.1$ and scale parameter $\beta=1.4$. It
is assumed that the system fails when a degradation process exceeds the
failure threshold $L=10$. The following costs are also assumed (m.u. monetary
units, t.u. time units )

$C_{p}=100\,m.u.,\quad C_{c}=200\,m.u.,\quad C_{I}=50\,m.u.,\quad
C_{d}=60\,m.u.\,per\,t.u..$ (27)

The search of the optimal maintenance policy corresponds to find the pair
$(T,M)$ that optimize the function $C(T,M)$ given by (25). That is, to find
$T_{opt}$ and $M_{opt}$ that fulfill

$C(T_{opt},M_{opt})=\inf\left\\{C(T,M),\ \ M\leq L\right\\}$

To  obtain Figure 7, 6,000 simulations were performed on a grid of size $10$
in the interval $(0,25)$ for $T$ and a grid of size $8$ in $(0,10)$ for $M$.
By visualization, the optimal values obtained for the time between inspections
and the preventive threshold are $T_{opt}=6.3333$ and $M_{opt}=6.1429$ with an
optimal expected cost of $35.3005$ monetary units per unit time.

Several sensitivity analysis were conducted to determine how the stochastic
intensity affects the proposed model and to analyze the robustness of the
solution when different parameters vary. Assuming that the degradation
processes start according to the intensity given by Eq. (26), assuming that
the system fails when the degradation level of a process exceeds the failure
threshold $L=10$ and with $\beta$ deterministic, Table 1 shows the influence
of the shape and scale parameters of the gamma process on the optimal expected
cost rate with a shaded grey scale when $\alpha$ and $\beta$ vary from $1$ to
$1.9$ increasing by $0.15$ at each step. Along with Table 1, Tables 2 and 3
show the optimal values of $T$ and $M$ for each combination of values of
$\alpha$ and $\beta$. Table 4 studies the influence of preventive and
corrective maintenance costs on the total optimal cost. To observe the
variation of the cost for preventive maintenance, 9 points have been
considered in the interval $[90,110]$, while the cost for corrective
maintenance has been inspected in the interval $[190,210]$, also with 9
points. This table has been obtained by fixing the optimal values for $T$ and
$M$ obtained previously.

Table 1 shows that $\alpha$ and $\beta$ have both effect on the optimal
expected cost. It is worth noting that more differences in the optimal
expected cost are found when the scale parameter varies. Similarly, as we can
see in Table 2, both parameters have effect on the optimal value for the time
between inspections $T$. However Table 3 does not show a clear influence of
the parameters in the case of the optimal value for the preventive threshold
$M$. The same can be said for Table 4 regarding to the influence of
maintenance costs in the model.

As it is shown in Table 1, the optimal expected cost decreases as $\alpha$ and
$\beta$ decreases. However, as it is shown in Table 2, the optimal time
between inspections decreases when $\alpha$ and $\beta$ increases. Higher
values for the scale and shape parameter of the gamma process mean a faster
deterioration hence inspections tend to be more frequent to reduce costs.

Figure 7: Expected cost rate versus $T$ and $M$. $\alpha$ $\beta$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 19.8048 | 21.4489 | 25.7375 | 27.8152 | 30.7940 | 30.7778 | 35.7045
1.15 | 22.4085 | 24.1170 | 28.1522 | 31.0347 | 32.3892 | 35.7728 | 39.1335
1.3 | 24.9069 | 27.3980 | 31.9920 | 33.3075 | 36.1599 | 39.5103 | 41.4842
1.45 | 27.6355 | 29.7365 | 34.2041 | 34.8008 | 37.3992 | 40.4628 | 46.5984
1.6 | 29.0795 | 30.4875 | 34.9262 | 38.3855 | 41.6051 | 46.3327 | 46.8082
1.75 | 30.4135 | 35.9786 | 37.9817 | 44.6625 | 45.7936 | 48.6405 | 52.9519
1.9 | 30.5336 | 36.1878 | 42.7939 | 45.5782 | 49.2110 | 52.4083 | 60.1414
Table 1: Sensitivity analysis on the values of $\alpha$ and $\beta$ for the optimal expected cost rate. $\alpha$ $\beta$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 8.2220 | 8.0979 | 6.7767 | 6.2052 | 6.0996 | 6.4665 | 4.5638
1.15 | 7.5315 | 7.1678 | 6.2315 | 5.7359 | 5.1743 | 5.2126 | 4.6818
1.3 | 7.4103 | 6.0552 | 6.5784 | 5.3184 | 5.0614 | 5.1942 | 4.8545
1.45 | 6.3486 | 5.9120 | 5.4764 | 5.2139 | 4.6110 | 4.5256 | 4.8004
1.6 | 5.8572 | 5.7846 | 5.2789 | 4.4831 | 4.1643 | 4.3130 | 3.9727
1.75 | 5.2285 | 4.9454 | 4.6661 | 4.3905 | 4.2302 | 3.9181 | 4.2534
1.9 | 5.6905 | 4.6343 | 4.3063 | 3.9870 | 3.9727 | 4.2653 | 4.7039
Table 2: Sensitivity analysis on the values $\alpha$ and $\beta$ for the optimal value of $T$. $\alpha$ $\beta$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 3.9426 | 2.1435 | 2.6667 | 1.6410 | 1.7250 | 2.7337 | 1.6285
1.15 | 1.2083 | 2.6261 | 2.1233 | 2.6294 | 2.6088 | 1.6989 | 2.9067
1.3 | 4.0309 | 1.3344 | 2.9083 | 2.5836 | 2.2619 | 2.0369 | 2.3512
1.45 | 1.1294 | 1.6428 | 1.3248 | 2.4431 | 1.8680 | 1.5133 | 1.8599
1.6 | 2.2714 | 1.8112 | 2.3257 | 1.3541 | 2.8780 | 1.4741 | 1.3602
1.75 | 1.9457 | 2.2737 | 1.1288 | 1.2540 | 2.0906 | 1.1675 | 1.9039
1.9 | 2.3188 | 2.8609 | 3.0805 | 1.2140 | 2.1427 | 1.6137 | 1.5923
Table 3: Sensitivity analysis on the values of $\alpha$ and $\beta$ for the optimal value of $M$. $C_{c}$ $C_{p}$ | 90 | 92.5 | 95 | 97.5 | 100 | 102.5 | 105 | 107.5 | 110
---|---|---|---|---|---|---|---|---|---
190 | 32.05 | 32.11 | 35.77 | 34.97 | 36.74 | 35.35 | 36.91 | 36.56 | 35.76
192.5 | 33.32 | 38.32 | 33.63 | 31.32 | 35.63 | 40.54 | 37.32 | 41.83 | 38.93
195 | 46.17 | 33.80 | 50.28 | 51.64 | 51.88 | 52.29 | 36.58 | 37.80 | 40.59
197.5 | 38.60 | 41.46 | 37.50 | 37.91 | 36.74 | 38.06 | 36.38 | 44.02 | 38.95
200 | 33.52 | 36.48 | 40.35 | 39.34 | 35.00 | 37.06 | 33.76 | 33.77 | 40.38
202.5 | 34.21 | 35.43 | 39.97 | 38.45 | 39.47 | 35.05 | 36.99 | 36.51 | 36.97
205 | 34.38 | 35.60 | 36.67 | 35.68 | 35.27 | 37.16 | 33.11 | 45.96 | 49.05
207.5 | 35.04 | 34.76 | 32.51 | 35.66 | 40.10 | 36.01 | 40.27 | 37.81 | 38.92
210 | 36.57 | 37.45 | 37.00 | 38.76 | 37.74 | 39.62 | 43.93 | 39.07 | 39.50
Table 4: Sensitivity analysis on the values of $C_{c}$ and $C_{p}$ for the
optimal expected cost rate.

###  7.2 Random effects model

This section assumes that the scale parameter $\beta$ is random. As before,
the degradation processes start according to the stochastic failure rate given
by Eq. (26) being $N(t)$ a Poisson process with parameter $\mu=2$. We assume
that the shape parameter of the gamma process is equal to $\alpha=1.1$ and the
scale parameter $\beta^{-1}$ follows a uniform distribution in the interval
$(1/1.4-\alpha^{*},1/1.4+\alpha^{*})$, with $\alpha^{*}=0.1$. As above, the
system fails when the level of degradation of a process exceeds $L=10$ and the
sequence of costs given by Eq. (27) is imposed. To analyse the optimal
maintenance strategy, a grid of size 10 in the interval $(0,25)$ has been
considered for $T$. Similarly, a grid of size 8 has been considered for $M$ in
$(0,10)$. For each combination of $M$ and $T$, 6,000 simulations were
performed. Figure 8 shows the results of the expected cost rate versus $T$ and
$M$. By inspection, the optimal values are obtained for $T_{opt}=6.3333$ and
$M_{opt}=4.8571$ with an optimal expected cost rate equals to $37.1962$
monetary units per  time unit.  Compared to the model without heterogeneity of
the previous subsection, both objective cost functions show the same pattern.
The case without heterogeneity shows a lower optimal expected cost rate
($C(T_{opt},M_{opt})=35.3005)$ monetary units per time unit than the case with
heterogeneity $C(T_{opt},M_{opt})=37.1962)$ monetary units per time unit (both
processes have the same expectation).

The sensitivity analysis of the expected cost rate and the optimal values for
$T$ and $M$ is summarized in Tables 5, 6 and 7, with $\alpha$ varying from $1$
to $1.9$ with increments of $0.15$ units and $\beta^{-1}$ following a uniform
distribution in $(b-\alpha^{*},b+\alpha^{*})$, being $\alpha^{*}=0.1$ and $b$
varying in the same manner as $\alpha$, from $1$ to $1.9$ also with increments
of $0.15$ units. Table 8 represents the variation in the expected cost rate
when the costs for preventive maintenance and corrective maintenance vary,
using the values of $T$ and $M$ obtained in the optimization of the model.

As in the model without heterogeneity, the optimal expected cost rate shown in
Table 5 also increases as $\alpha$ increases. That is, more deterioration
implies higher costs. In this particular example, the random choice of the
parameter $\beta$ has the same effect as $\alpha$ on the resulting expected
cost rate: it increases when $b$ increases. Similarly, Table 6 shows that the
optimal value of $T$ decreases as $\alpha$ increases. An increase in the
deterioration implies more frequent inspections. Finally, no trend seems to be
found for the optimal values of the preventive threshold $M$ or in the
variation of maintenance costs $C_{c}$ and $C_{p}$ shown in Tables 7 and 8,
respectively.

Figure 8: Expected cost rate versus $T$ and $M$ in a model with random effects. $\alpha$ $b$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 26.1538 | 28.3350 | 31.8173 | 32.9500 | 38.0842 | 44.5595 | 48.0884
1.15 | 27.2921 | 30.9587 | 34.4670 | 38.8324 | 45.0903 | 46.3951 | 54.5271
1.3 | 37.6802 | 37.1068 | 42.1956 | 44.4267 | 49.5072 | 56.1529 | 58.2646
1.45 | 32.5948 | 39.7887 | 44.3998 | 49.9204 | 53.8380 | 56.2251 | 61.2970
1.6 | 37.4109 | 42.5415 | 48.5337 | 53.8215 | 55.9734 | 60.1671 | 67.9211
1.75 | 39.4109 | 46.8386 | 56.7378 | 58.7266 | 60.9805 | 61.7760 | 64.8222
1.9 | 46.9104 | 51.9782 | 57.3185 | 60.1173 | 63.9612 | 66.7328 | 66.8821
Table 5: Sensitivity analysis on the values of $\alpha$ and $b$ for the optimal expected cost rate in a model with random effects. $\alpha$ $b$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 8.3615 | 9.3511 | 10.8538 | 11.6696 | 13.6147 | 13.8070 | 14.9346
1.15 | 7.4953 | 8.3954 | 9.3891 | 10.4470 | 11.4916 | 12.3747 | 13.6643
1.3 | 6.7435 | 8.0170 | 8.8949 | 9.2960 | 10.3069 | 11.4964 | 12.1105
1.45 | 6.5243 | 6.9606 | 7.9728 | 9.0949 | 9.4303 | 10.1668 | 11.6343
1.6 | 5.8121 | 6.7253 | 8.0823 | 8.6948 | 9.0396 | 9.4426 | 10.3791
1.75 | 5.7501 | 6.0871 | 6.9204 | 7.4151 | 8.1554 | 8.7310 | 9.5734
1.9 | 5.4413 | 5.8472 | 6.4068 | 6.8895 | 8.0662 | 8.1205 | 8.8067
Table 6: Sensitivity analysis on the values of $\alpha$ and $b$ for the optimal value of $T$ in a model with random effects. $\alpha$ $b$ | 1 | 1.15 | 1.3 | 1.45 | 1.6 | 1.75 | 1.9
---|---|---|---|---|---|---|---
1 | 2.7979 | 1.1639 | 2.2698 | 3.6165 | 1.6126 | 3.4845 | 3.3828
1.15 | 1.3415 | 5.5418 | 1.6695 | 3.1643 | 1.8576 | 2.7601 | 3.7015
1.3 | 1.4256 | 1.1350 | 1.6613 | 3.2311 | 3.1447 | 1.9884 | 2.3213
1.45 | 2.9804 | 1.7815 | 1.5055 | 1.4391 | 2.4255 | 1.5748 | 5.2194
1.6 | 1.7922 | 2.3719 | 3.4496 | 2.5720 | 3.2729 | 3.6388 | 1.6719
1.75 | 3.5678 | 2.5439 | 3.6870 | 1.2555 | 3.2687 | 2.2337 | 3.3523
1.9 | 1.4376 | 3.1995 | 1.1186 | 2.0668 | 1.3907 | 4.2251 | 2.5016
Table 7: Sensitivity analysis on the values of $\alpha$ and $b$ for the optimal value of $M$ in a model with random effects. $C_{c}$ $C_{p}$ | 90 | 92.5 | 95 | 97.5 | 100 | 102.5 | 105 | 107.5 | 110
---|---|---|---|---|---|---|---|---|---
190 | 31.15 | 36.48 | 33.77 | 36.77 | 33.09 | 33.18 | 43.31 | 39.84 | 33.40
192.5 | 37.24 | 35.31 | 32.65 | 33.58 | 37.47 | 33.56 | 34.54 | 32.91 | 31.88
195 | 32.96 | 33.19 | 34.71 | 31.13 | 32.93 | 31.62 | 34.08 | 34.78 | 33.02
197.5 | 33.65 | 45.19 | 34.26 | 36.21 | 34.43 | 43.83 | 36.55 | 42.48 | 34.34
200 | 33.51 | 43.01 | 39.67 | 31.93 | 39.48 | 36.13 | 34.60 | 33.81 | 37.08
202.5 | 34.02 | 36.24 | 33.92 | 35.40 | 36.48 | 33.11 | 35.94 | 34.56 | 34.61
205 | 35.99 | 33.67 | 44.57 | 38.59 | 35.89 | 39.22 | 33.87 | 41.78 | 38.03
207.5 | 34.59 | 30.98 | 33.59 | 35.03 | 36.18 | 43.77 | 50.11 | 36.11 | 33.27
210 | 39.02 | 34.31 | 37.23 | 34.18 | 34.70 | 35.04 | 36.52 | 39.06 | 39.01
Table 8: Sensitivity analysis on the values of $C_{c}$ and $C_{p}$ for the
optimal expected cost rate in a model with random effects.

##  8 Conclusions, future extensions and limitations of this work

A  system subject to multiple degradation processes is analyzed in this paper.
Degradation processes start according to a shot noise Cox process and grow
according to homogeneous gamma processes Under this framework, the combined
model of initiation and growth is modelled as a Cox process. Using properties
of the Cox processes, the distribution to the system lifetime is obtained. It
is shown that, in absence of maintenance, the system lifetime  distribution is
increasing failure rate hence a preventive maintenance policy is worth
implementing to improve the system reliability. The analysis is also completed
with degradation process-specific heterogeneity considering a random effects
model  using a uniform distribution. A classical maintenance strategy is
implemented for this system. As before, the probability expressions of the
objective cost function are closed forms that allow intuitive interpretations.
It is due to the mathematical tractability of the shot noise Cox process
induced by a Poisson process.

Although in this paper starting points of the degradation processes follow a
shot noise Cox process, the result can be extended considering a different Cox
process such as a Weibull renewal process. Or, even, a non-Cox distribution
such as a Hawkes process. Another crucial assumption in this paper is that the
degradation processes evolve independently and according to the same
degradation pattern. For future research, different degradation patters for
the degradation processes and dependence between the processes could be
considered.

With respect to the limitations of the paper, the most important limitation
refers to the lack of real data that support the proposed maintenance model.

## Acknowledgements

This research was supported by Ministerio de Ciencia e Innovación, Spain
(Project PGC2018-094964-B-I00).

## References

## References

  * [1] Alaswad S., Xiang Y. A review on condition-based maintenance optimization models for stochastically deteriorating system, Reliability Engineering $\&$ System Safety, 2017; 157, 54–63.
  * [2] Bautista L., Castro I.T., Landesa L. Maintenance cost assessment for heterogeneous multi-component systems incorporating perfect inspections and waiting time to maintenance. To appear in Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability https://doi.org/10.1177/1748006X211038804
  * [3] Bertoin J. Lévy processes. Cambridge University Press, 1996.
  * [4] Bordes L., Mercier S., Remy E., Dautreme E. Partially observed competing degradation processes: modeling and inference, Applied Stochastic Models in Business and Industry, 2016; 32(5), 555-739.
  * [5] Caballé N., Castro I.T., Pérez C.J., Lanza-Gutiérrez J. A condition-based maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes, Reliability Engineering $\&$ System Safety, 2015; 134, 89-109.
  * [6] Castro I.T., Caballé N., Pérez C.J. A condition-based maintenance for a system subject to multiple degradation processes and external shocks, International Journal of Systems Science, 2015; 46(9), 1692-1704.
  * [7] Castro I.T., Mercier S. Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, 2016; 230(4), 364–377.
  * [8] Castro I.T., Landesa L. A dependent complex degrading system with non-periodic inspection times, Computers $\&$ Industrial Engineering, 2019; 133, 241–252.
  * [9] Cha J.H., Finkelstein M. Stochastic survival models with events triggered by external shocks, Probability in the Engineering and Informational Sciences, 2012; 26, 183-195.
  * [10] Cha J.H., Finkelstein M. On Some Shock Models with Poisson and Generalized Poisson Shock Processes. ICSA Book Series in Statistics, 2017; 67–79.
  * [11] Cha J.H., Finkelstein M. On information-based residual lifetime in survival models with delayed failures, Statistic $\&$ Probability Letters, 2018; 137, 209-216.
  * [12] Cha J.H., Finkelstein M. On a new shot noise process and the induced survival model, Methodology and Computing in Applied Probability, 2018; 20, 897-917.
  * [13] Cha J.H., Finkelstein M, Levitin G. Bivariate preventive maintenance of systems with lifetimes dependent on a random shock process, European Journal of Operational Research, 2018; 266, 122-143.
  * [14] Chen, N., Ye Z.S., Xiang Y., Zhang L. Condition based maintenance using the inverse Gaussian degradation model, European Journal of Operational Research, 2015; 243(1), 190-199.
  * [15] Chen, J. M., Hawkes, A. G., Scalas, E., Trinh, M. Performance of information criteria for selection of Hawkes process models of financial data, Quantitative Finance, 2018; 18(2), 225–235.
  * [16] Crk V. Reliability assessment from degradation data. 2000 Proceedings Annual Reliability and Maintainability Symposium, 2000; 155-161.
  * [17] Finkelstein M., Cha J.H., Levitin G. A hybrid preventive maintenance model for systems with partially observable degradation, IMA Journal of Management Mathematics, 2020; 31(3), 345-365.
  * [18] Guida M., Calabria R., Pulcini G. Bayes inference for a non-homogeneous Poisson process with power intensity law (reliability), IEEE Transactions on Reliability, 1989; 38(5), 603–609.
  * [19] Huang Y., Vasan A.S.S., Doraiswami R., Osterman M. MEMS Reliability Review, IEEE Transactions on Device and Materials Reliability, 2012; 12(2), 482–493.
  * [20] Huynh K.T., Grall A., Bérenguer C. Assessment of diagnostic and prognostic condition indices for efficient and robust maintenance decision-making of systems subject to stress corrosion cracking, Reliability Engineering $\&$ System Safety, 2017; 159, 237–254.
  * [21] Huynh K.T. Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems, European Journal of Operational Research, 2020; 280(1), 152-163.
  * [22] Ishiara S., Goshima T., Nomura K., Yoshimot T. Crack propagation behavior of cermets and cemented carbides under repeated thermal shocks by the improved quench test, Journal of Materials Science, 1999; 34, 629–636.
  * [23] Jia G., Gardoni P. State-dependent stochastic models: A general stochastic framework for modeling deteriorating engineering systems considering multiple deterioration processes and their interactions, Structural Safety, 2018; 72, 99–110.
  * [24] Kuniewski S.P., van der Weide J.A.M., van Noortwijk J.M. Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection, Reliability Enginering $\&$ System Safety, 2009; 94, 1480–1490.
  * [25] Lawless J., Crowder M. Covariates and random effects in a gamma process model with application to degradation and failure, Lifetime Data Analysis, 2004; 10, 213-227.
  * [26] Lemoine A.J., Wenocur M.L. On failure modeling, Naval Research Logistics Quarterly, 1985; 32(3), 497-508.
  * [27] Lemoine A.J., Wenocur M.L. A note on shot-noise and reliability modeling, Operations Research, 1986; 34(2), 320-323.
  * [28] Liu B., Zhao X., Liu G., Liu Y. Life cycle cost analysis considering multiple dependent degradation processes and environmental influence, Reliability Engineering $\&$ System Safety, 2020; 197, 106784.
  * [29] Möller J. Properties of spatial Cox process models, Journal of Statistical Research of Iran, 2005; 2(1), 1–18.
  * [30] Pulcini G. Modeling the mileage accumulation process with random effects, Communications in Statistics-Theory and Methods, 2013; 42(15), 2661-2683.
  * [31] Qiu Q., Cui L., Yang L. Maintenance policies for energy systems subject to complex failure processes and power purchasing agreement, Computers $\&$ Industrial Engineering, 2018; 119, 193-203.
  * [32] Rodríguez-Picón L.A., Rodríguez-Picón A.P., Méndez-González L.C., Rodríguez-Borbón M.I., Alvarado-Iniesta A. Degradation modeling based on gamma process models with random effects, Communications in Statistics - Simulation and Computation, 2018; 47(6), 1796-1810.
  * [33] Ross S.M. Introduction to Probability Models, 2014, Elsevier, Rome.
  * [34] Serfozo R.F. Point processes. Handbooks in Operations Research and Management Science. Vol 2: Stochastic Models (Heyman D.P. $\&$ Sobel M.J. eds) (1990) 1–93. North Holland, Amsterdam.
  * [35] Tee K.F., Pesinis K. Reliability prediction for corroding natural gas pipelines, Tunnelling and Underground Space Technology, 2017; 65, 91–85.
  * [36] Tijms H.C. A First Course in Stochastic Models. 2013, John Wiley $\&$ Sons Ltd, West Sussex (England).
  * [37] Tsai C.C., Tseng S.T., Balakrishnan N. Optimal design for degradation tests based on gamma processes with random effects, EEE Transactions on Reliability, 2012; 61(2), 604-613.
  * [38] Wang Y., Pham H. Modeling the dependent competing risks with multiple degradation processes and random shocks using time-varying copulas, IEEEE Transactions on Reliability, 2012; 61, 13-22.
  * [39] Wu S., Castro I.T. Maintenance policy for a system with a weighted linear combination of degradation processes, European Journal of Operational Research, 280(1), 2020; 124-133.
  * [40] Xiang Y., Zhu Z., Coit D.W., Feng Q. Condition-based maintenance under performance-based contracting, Computers $\&$ Industrial Engineering, 2017; 111, 391-402.