# Assessment of the maintenance cost and analysis of availability measures in
a finite life cycle for a system subject to competing failures

N. C. Caballé  I. T. Castro<EMAIL_ADDRESS>Department of Mathematics,
University of Extremadura, Spain

###### Abstract

This paper deals with the assessment of the performance of a system under a
finite planning horizon. The system is subject to two dependent causes of
failure: internal degradation and sudden shocks. We assume that internal
degradation follows a gamma process. When the deterioration level of the
degradation process exceeds a predetermined value, a degradation failure
occurs. Sudden shocks arrive at the system following a doubly stochastic
Poisson process (DSPP). A sudden shock provokes the total breakdown of the
system. A condition-based maintenance (CBM) with periodic inspection times is
developed. To evaluate the maintenance cost, recursive methods combining
numerical integration and Monte Carlo simulation are developed to evalute the
expected cost rate and its standard deviation. Also, recursive methods to
calculate some transient measures of the system are given. Numerical examples
are provided to illustrate the analytical results.

###### keywords:

Availability , condition based maintenance , degradation threshold shock model
, interval reliability , gamma process , reliability.

††journal: ArXiv

## Acronyms

CBM: condition-based maintenance.
CM: corrective maintenance.
DSPP: doubly stochastic Poisson process.
DTS: degradation-threshold-shock.
$d.u.$: degradation units.
HPP: homogeneous Poisson process.
$m.u.$: monetary units.
NHPP: non-homogeneous Poisson process.
PM: preventive maintenance.
$t.u.$: time units.

## Notation

$A_{T}^{M}(t)$: availability of the system at time $t$ with time between
inspections $T$ and preventive threshold $M$.

$\tilde{A}_{T}^{M}(t)$: estimation of $A_{T}^{M}(t)$.

$\alpha$, $\beta$: shape and scale parameters of the gamma process.

$\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}$: gamma process parameters modified
under the $i$-th and $j$-th positions for the vector $\mathbf{v}$.

$C_{c}$: corrective maintenance cost.

$C_{d}$: downtime cost.

$C_{I}$: inspection cost.

$C_{p}$: preventive maintenance cost.

$C(t_{1},t_{2})$: maintenance cost in $(t_{1},t_{2}]$.

$C(t_{f})$: maintenance cost in the life cycle.

$C^{\infty}(T,M)$: expected cost rate for a time between inspections $T$ and
preventive threshold $M$.

$D_{j}$: length of the $j$-th renewal cycle.

$E\left[C_{T}^{M}(t)\right]$: expected transient cost at time $t$ with time
between inspections $T$ and preventive threshold $M$.

$\tilde{E}\left[C_{T}^{M}(t)\right]$: estimation of
$E\left[C_{T}^{M}(t)\right]$.

$E\left[C_{T}^{M}(t)^{2}\right]$: expected square cost at time $t$ with time
between inspections $T$ and preventive threshold $M$.

$\tilde{E}\left[C_{T}^{M}(t)^{2}\right]$: estimation of
$E\left[C_{T}^{M}(t)^{2}\right]$.

$\tilde{E}^{*}\left[C^{M}_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}(t_{f})\right]$:
minimal expected transient cost varying gamma process parameters.

$\tilde{E}^{*}\left[C^{M}_{T,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(t_{f})\right]$:
minimal expected transient cost varying sudden shock process parameters.

$E\left[W_{T}^{M}(t_{1},t_{2})\right]$: expected downtime in $(t_{1},t_{2}]$
for a time between inspections $T$ and preventive threshold $M$.

$\tilde{E}\left[W_{T}^{M}(t_{1},t_{2})\right]$: estimation of
$E\left[W_{T}^{M}(t_{1},t_{2})\right]$.

$f_{\alpha t,\beta}$: density function of the gamma process.

$F_{\sigma_{z}}$: distribution function of $\sigma_{z}$.

$\bar{F}_{\sigma_{z_{2}}-\sigma_{z_{1}}}$: survival function of
$\sigma_{z_{2}}-\sigma_{z_{1}}$.

$I(v,t)$: survival function of $Y$ conditioned to $\sigma_{M_{s}}=v$.

$IR_{T}^{M}(t_{1},t_{2})$: interval reliability in $(t_{1},t_{2}]$ with time
between inspections $T$ and preventive threshold $M$.

$L$: breakdown threshold.

$\lambda(t)$: intensity of the sudden shock process.

$\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}$: sudden shock process parameters
modified under the $i$-th and $j$-th positions of the vector $\mathbf{v}$.

$M$: preventive threshold.

$M_{s}$: threshold from which the system is more prone to sudden shocks.

$N(t)$: number of renewals up to $t$.

$N_{s}(t_{1},t_{2})$: number of sudden shocks in $(t_{1},t_{2}]$.

$O(t)$: deterioration level of the maintained system at time $t$.

$P_{R_{1}}^{M}(kT)$: probability of a maintenance action at $k$-th inspection
for a time between inspections $T$ and preventive threshold $M$.

$\tilde{P}_{R_{1}}^{M}(kT)$: estimation of $P_{R_{1}}^{M}(kT)$.

$P_{R_{1,p}}^{M}(kT)$, $P_{R_{1,c}}^{M}(kT)$: probability of a preventive and
corrective maintenance action at $k$-th inspection for a time between
inspections $T$ and preventive threshold $M$, respectively.

$\tilde{P}_{R_{1,p}}^{M}(kT)$, $\tilde{P}_{R_{1,c}}^{M}(kT)$: estimation of
$P_{R_{1,p}}^{M}(kT)$ and $P_{R_{1,c}}^{M}(kT)$.

$R_{j}$: chronological time of the $j$-th renewal cycle.

$R_{T}^{M}(t)$: reliability of the system at time $t$ with time between
inspections $T$ and preventive threshold $M$.

$\tilde{R}_{T}^{M}(t)$: estimation of $R_{T}^{M}(t)$.

$S_{T}^{M}(t)$: standard deviation of the expected cost at time $t$ with time
between inspections $T$ and preventive threshold $M$.

$\sigma_{z}$: time to reach the deterioration level $z$.

$T$: time between inspections.

$t_{f}$: length of the life cycle.

$\mathbf{v}$: variation vector for model parameters.

$V^{M}_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}(t_{f})$: relative variation
percentage for gamma process parameters.

$V^{M}_{T,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(t_{f})$ relative
variation percentage for sudden shock process parameters.

$X(t)$: underlying degradation process.

$Y$: time to a sudden shock.

##  1 Introduction

A fundamental aim in the industry field is to ensure the reliability of the
systems. It is well-known that some systems suffer a physical degradation
process which precedes the failure. This degradation process may involve
chemical and physical changes in the system complicating its maintenance. The
theory of stochastic processes provides an analytical framework for modelling
the impact of the uncertain and time-dependent degradation processes.

The gamma process is a stochastic cumulative process considered as one of the
most appropriated processes for modelling the damage involved by the
cumulative deterioration of systems and structures [1]. It is characterised by
independent and non-negative gamma increments with identical scale parameters.
The gamma process was first applied by Moran [2] to model water flow into a
dam. Later, Abdel-Hammed [3] proposed the gamma process as a specific model
for deterioration occurring randomly in time. From then on, gamma process has
been widely used in the reliability field. The survey by Van Noortwijk [1]
provides many examples of the use of the gamma process in engineering.

However, some systems are not only subject to internal degradation but are
also exposed to sudden shocks which can cause its failure. For example, an
ammeter degrades over time and also receive shocks (such as ligthtning) that
could provoke the failure of the ammeter [4]. The Light-Emitting Diode (LED)
lamps are subject to underlying degradation processes and shocks (over-voltage
and over-heating) [5]. Even the human body can be regarded as a competing risk
system where some health markers can degrade due to the age and some shocks
(catastrophic or not) can happen (as a heart attack). Up to our knowledge,
Lemoine and Wenocur [6] were the first to combine both causes of failure,
proposing the Degradation-Threshold-Shock (DTS) models. Singpurwalla [7]
describes diverse failure models for DTS models. Since two competing failure
causes are analyzed, the analysis can be performed considering independent
causes of failure [8] or dependent [9] [10],[11] .

Maintenance strategies regulate the different maintenance tasks which must be
performed on the system. Establishing a good maintenance task, the correct
functioning of the system is ensured and the maintenance cost can be
optimised. Condition-based maintenance (CBM) is one of the most popular
techniques used for degrading systems. CBM is a maintenance program that
recommends to perform maintenance actions based on the information collected
through a condition monitoring process using certain types of sensors or other
appropriate indicators [12]. The implementation of CBM programs for DTS models
is not new (see [13] and [14] among others)

Many maintenance designs are planned based on an infinite operating horizon.
It means that, after any replacement, the system is renewed by a new one with
the same characteristics and the same process is assumed to be repeated
indefinitely. Characteristics of these systems, such as the degradation level
or the age, are often selected as criteria to optimise the long-run cost rate.
Due to renewal properties, this long-run cost rate is equal to the expected
cost in a renewal cycle divided by the length of the renewal cycle (see e.g.
[15, 16, 17, 18]). However, most systems actually have a finite operating life
cycle since the system cannot always be replaced by a new one with the same
characteristics as the previous one an infinite number of times. For instance,
in military applications, a missile launching system is only required to be
functioning within the designated mission time [19]. Hence, the use of the
asymptotic approach is questionable and the maintenance cost should be
analysed under a transient approach.

Although the transient approach is more realistic than the asymptotic
approach, it is less used due to the analytical and computational difficulty
of treatment that it involves. However, some works can be found in the
literature where authors analyse maintenance strategies under a transient
approach for systems subject to competing risks. For example, Taghipour et al.
[20] proposed a model to find the optimal interval periodic inspection
interval on a finite life cycle for a system subject to different types of
failure.

This paper expands the works by Cheng et al. [21] and Pandey et al. [22]
considering the time as a continuous variable and by adding a new component of
risk (sudden shocks), whose arrival depends on the degradation process of the
system. The framework exposed in [25] inspires the setting of this problem. In
this paper, we assume that the system is degraded following a gamma process
and sudden shocks arrive at the system following a doubly stochastic Poisson
process (DSPP) whose intensity function depends on the degradation process. A
CBM with periodical inspection times is developed using the expected cost rate
in the life cycle as objective cost function. The evaluation of the
maintenance cost in the life cycle of the system is performed using recursive
methods. Furthermore, the results obtained using recursive methods are
compared to the results obtained based on strictly Monte Carlo simulation.
Further comparisons of the maintenance cost are performed considering an
infinite life cycle. The robustness of the gamma process parameters and the
shock process is also analysed.

In many application fields, there is an increasing interest in evaluating the
performance of maintained systems. Reliability and availability are two
important performance measures in the traditional reliability field. But some
situations are not covered by these indexes and new performance measures, such
as the interval reliability, has been developed [23, 24]. Along with the
maintenance assessment, in this paper, some performance measures are also
evaluated in the life cycle of the system. The evaluation of these performance
measures in the life cycle of the system is performed using recursive methods.
Furthermore, the results obtained using recursive methods are compared to the
results obtained based on strictly Monte Carlo simulation.

In short, the main contributions of this paper are:

  1. 1.

Development of a recursive method to obtain the expected cost and its standard
deviation in the life cycle of the system.

  2. 2.

Comparing the results obtained using this recursive method to the results
obtained by using strictly Monte Carlo simulation.

  3. 3.

Comparing the asymptotic expected cost rate and the expected cost in the life
cycle of the system.

  4. 4.

Analysis of the robustness of some parameters that describe the functioning of
the system.

  5. 5.

Assessment of the the availability, reliability, and interval reliability in
the life cycle of the system using a recursive method and comparison of the
results obtained by using strictly Monte Carlo simulation.

In order to develop these contributions, this paper is structured as follows.
In Section 2 the general framework of the model is described. Expected
transient cost and its standard deviation associated are analysed in Section
3. Performance measures of the system are exposed in Section 4. Numerical
examples are given in Section 5. Conclusions and further possible extensions
of this paper are provided in Section 6.

##  2 Framework of the problem

A system subject to two dependent competing causes of failure, degradation and
sudden shocks, is considered in this paper. The general assumptions of this
model are:

  1. 1.

The system starts working at time $t=0$. This system is subject to an internal
degradation process which evolves according to a homogeneous gamma process
with parameters $\alpha$ and $\beta,$ where $\alpha,\beta>0$. Let $X(t)$ be
the deterioration level of the system at time $t$. Thus, for two time instants
$s$ and $t$, with $s<t$, the density function of the increment deterioration
level $X(t)-X(s)$ is given by

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

where $\Gamma(\cdot)$ denotes the gamma function defined as

$\Gamma(\alpha)=\displaystyle\int_{0}^{\infty}u^{\alpha-1}e^{-u}\leavevmode\nobreak\
du.$ (2)

The system fails due to degradation when the deterioration level exceeds a
fixed threshold $L$, called the breakdown threshold.

  2. 2.

The system not only fails due to internal degradation, but also it is subject
to sudden shocks which cause its failure. Sudden shocks arrive at the system
according to a process $\left\\{N_{s}(t),t\geq 0\right\\}$. This process shows
the dependence between degradation and shocks. Following the spirit showed in
[25], we assume that $\left\\{N_{s}(t),t\geq 0\right\\}$ is a DSPP with
intensity $\lambda(t,X(t))$ given by

$\lambda\left(t,X(t)\right)=\lambda_{1}(t)\mathbf{1}_{\left\\{X(t)\leq
M_{s}\right\\}}+\lambda_{2}(t)\mathbf{1}_{\left\\{X(t)>M_{s}\right\\}},$ (3)

for $t\geq 0$, where $\lambda_{1}$ and $\lambda_{2}$ denote two failure rate
functions which verify $\lambda_{1}(t)\leq\lambda_{2}(t)$, for all $t\geq 0$
and where $\mathbf{1}_{\\{\cdot\\}}$ denote the indicator function which
equals $1$ if the argument is true and $0$ otherwise. The arrival of a sudden
shock provokes the system failure.

  3. 3.

The system is inspected each $T\ (T>0)$ time units ($t.u.$). In these
instants, it is checked if the system is working or is down. If the system is
down, a corrective maintenance (CM) is performed and the system is replaced by
a new one. A CM event simulation is shown in Fig. 1.

Fig. 1: A corrective maintenance event.

Let $M$ be the deterioration level from which the system is considered as too
worn ($M<L$). If the system is still working and the deterioration level of
the system exceeds the preventive threshold $M$, a preventive maintenance (PM)
is performed and the system is replaced by a new one. Otherwise, no
maintenance action is performed. A PM event simulation is shown in Fig. 2.

Fig. 2: A preventive maintenance event.

We assume that the time required to perform a maintenance action is
negligible.

  4. 4.

All maintenance actions imply a cost. A CM and a PM have associated a cost of
$C_{c}$ and $C_{p}$ monetary units ($m.u.$), respectively, and each inspection
implies a cost of $C_{I}$ $m.u.$ In addition, if the system fails, the system
is down until the next inspection. Each time unit that the system is down, a
cost of $C_{d}$ $m.u./t.u.$ is incurred. We assume $C_{c}>C_{p}>C_{I}$.

  5. 5.

Let $(0,t_{f}]$ be the finite operating life cycle of the system. It means
that, if the calendar time exceeds $t_{f}$, the system can no longer be
replaced by a new one with the same charasteristics.

Let $\sigma_{z}$ be the random variable describing the time to reach a certain
degradation level $z$. The distribution function of $\sigma_{z}$, denoted as
$F_{\sigma_{z}}$, is given by

$\begin{array}[]{l}\begin{aligned} F_{\sigma_{z}}(t)=&P[X(t)\geq z]\\\
=&\displaystyle\int_{z}^{\infty}f_{\alpha t,\beta}(x)\leavevmode\nobreak\
dx=\frac{\Gamma(\alpha t,z\beta)}{\Gamma(\alpha t)},\end{aligned}\end{array}$
(4)

for $t\geq 0$ where $f_{\alpha t,\beta}(x)$ and $\Gamma(\alpha t)$ are given
by (1) and (2), respectively, and

$\Gamma(\alpha,x)=\displaystyle\int_{x}^{\infty}u^{\alpha-1}e^{-u}\leavevmode\nobreak\
du,$

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

For $z_{1}\leq z_{2}$, the survival function of
$\sigma_{z_{2}}-\sigma_{z_{1}}$ is given by

$\begin{array}[]{l}\begin{aligned}
\bar{F}_{\sigma_{z_{2}}-\sigma_{z_{1}}}(t)=&P\left[\sigma_{z_{2}}-\sigma_{z_{1}}\geq
t\right]\\\
=&\displaystyle\int_{x=0}^{\infty}\int_{y=z_{1}}^{\infty}f_{\sigma_{z_{1}},X(\sigma_{z_{1}})}(x,y)\\\
&F_{\alpha t,\beta}(z_{2}-y)\leavevmode\nobreak\ dy\leavevmode\nobreak\
dx,\end{aligned}\end{array}$ (5)

where $F_{\alpha t,\beta}$ denotes the distribution function of $f_{\alpha
t,\beta}$ and $f_{\sigma_{z_{1}},X(\sigma_{z_{1}})}$ denotes the joint density
function of $\left(\sigma_{z_{1}},X(\sigma_{z_{1}})\right)$ provided by
Bertoin [26] as

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

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

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

From Assumption 2, the shock process follows a doubly stochastic Poisson
process where the intensity of the shocks depends on time $t$ and on the
degradation level $X(t)$. It means that, given a path $x$ de $X(t)$, the
process $\left\\{N_{s}(t),t\geq 0\right\\}$ is a non homogeneous Poisson
process with intensity $\lambda(t,x)$. In absence of maintenance, the time to
a failure system is defined as the minimum $D=\min(\sigma_{L},Y)$ where

$Y=\inf\left(t\geq 0,N_{s}(t)=1\right),$

with survival function

$P(D>t)=\mathbb{E}\left[\mathbf{1}_{\left\\{D>t\right\\}}exp\left(-\int_{0}^{t}\lambda(s,X(s)ds)\right)\right]$

Let $I(v,t)$ be the survival function of $Y$ for $t\geq v$, conditioned to
$\sigma_{M_{s}}=v$. That is

$\begin{array}[]{l}\begin{aligned}
I(v,t)=&P\left[Y>t|_{\sigma_{M_{s}}=v}\right]\\\
=&\exp\left\\{-\displaystyle\int_{0}^{t}\lambda\left(z\right)\leavevmode\nobreak\
dz\right\\}=\frac{\bar{F}_{1}(v)}{\bar{F}_{1}(0)}\frac{\bar{F}_{2}(t)}{\bar{F}_{2}(v)},\end{aligned}\end{array}$
(6)

where

$\bar{F}_{j}(t)=\exp\left\\{-\displaystyle\int_{0}^{t}\lambda_{j}(u)du\right\\},\quad
j=1,2,$ (7)

with density function $f_{j}(t)$, for $j=1,2$.

##  3 Expected transient cost analysis

A goal in industry is to find the maintenance strategy that minimises an
objective cost function. One of the most used function in the literature as
objective cost function is the asymptotic cost rate [22] (asymptotic cost per
time unit) that has a simple expression if the functioning of the system can
be modelled as a renewal process.

Since in this paper the system is repaired after each preventive or corrective
maintenance, let $D_{1},D_{2},\ldots$ be the time to the successive renewals
of the system. Let $C^{\infty}(T,M)$ be the asymptotic cost rate with a time
between inspections $T$ and preventive threshold $M$. Based on the “Renewal
Theorem”, $C^{\infty}(T,M)$ is equal to the expected cost in a renewal cycle
divided by the length of the renewal cycle. That is

$C^{\infty}(T,M)=\lim\limits_{t\rightarrow\infty}\frac{C(t)}{t}=\frac{E[C_{1}]}{E[D_{1}]},$

where $C(t)$ denotes the maintenance cost at time $t$, and $C_{1}$ and $D_{1}$
the cost and the length of a renewal cycle, respectively. In this paper,
$C^{\infty}(T,M)$ is given by

$\begin{array}[]{l}\begin{aligned}
&\frac{\displaystyle\sum_{k=1}^{\infty}\Big{[}C_{c}P_{R_{1},c}^{M}(kT)+C_{p}P_{R_{1},p}^{M}(kT)+C_{I}(k-1)P_{R_{1}}^{M}(kT)\Big{]}}{\displaystyle\sum_{k=1}^{\infty}kT\leavevmode\nobreak\
P_{R_{1}}^{M}(kT)}\\\
+&\frac{\displaystyle\sum_{k=1}^{\infty}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]}{\displaystyle\sum_{k=1}^{\infty}kT\leavevmode\nobreak\
P_{R_{1}}^{M}(kT)},\\\ \end{aligned}\end{array}$ (8)

where $P_{R_{1}}^{M}(kT)$ denotes the probability of the first maintenance
action at time $kT$ for $k=1,2,\ldots$ given by

$P_{R_{1}}^{M}(kT)=P_{R_{1,p}}^{M}(kT)+P_{R_{1,c}}^{M}(kT),$ (9)

$P_{R_{1,p}}^{M}(kT)$ and $P_{R_{1,c}}^{M}(kT)$ denote the probability of the
first preventive and corrective maintenance action at time $kT$ for
$k=1,2,\ldots$, respectively, and $E\left[W_{T}^{M}((k-1)T,kT)\right]$ the
expected downtime in $((k-1)T,kT]$. The analytical expressions for these
quantities were provided by Huynh et al. [25].

Since in this paper the objective is to evaluate the maintenance cost in the
life cycle of the system, the expected cost rate in this life cycle is used as
objective cost function. Let $R_{j}$ be the chronological time of the $j$-th
renewal cycle, for $j=0,1,2,\ldots,N(t_{f})+1$, being $N(t_{f})$ the number of
complete renewals in the life cycle of the system. Thus, $R_{j}$ is given by

$R_{j}=\displaystyle\sum_{n=1}^{j}D_{n},$

for $R_{0}=0$, where $D_{n}$ denotes the length of the $n$-th renewal cycle,
with $n=1,2,\ldots,N(t_{f})+1$. Hence, the length of the $n$-th renewal cycle
$D_{n}$ is given by

$D_{n}=\left\\{\begin{array}[]{ll}R_{n}-R_{n-1},&\text{if}\quad
n=1,2,\ldots,N(t_{f})\\\ t_{f}-R_{N(t_{f})},&\text{if}\quad
n=N(t_{f})+1\end{array}\right..$

Fig. 3 shows a process realisation.

Fig. 3: A realization of a renewal cycle sequence.

The total cost in the finite life cycle is the sum of the incurred costs in
the different $N(t_{f})$ renewal cycles and the incurred cost in
$(R_{N(t_{f})},t_{f}]$. That is

$C(t_{f})=\displaystyle\sum_{j=1}^{N(t_{f})}C(R_{j-1},R_{j})+C(R_{N(t_{f})},t_{f}),$

where $C(t_{1},t_{2})$ denotes the cost in the interval $(t_{1},t_{2}]$, and
$C(0,t_{f})$ is simplified as $C(t_{f})$.

Let $E\left[C_{T}^{M}(t)\right]$ be the expected transient cost at time $t$
with a time between inspections $T$ and a preventive threshold $M$. The next
result provides the Markov renewal equation that fulfils
$E\left[C_{T}^{M}(t)\right]$.

######  Theorem 1.

For $t<T$, the expected transient cost, $E\left[C_{T}^{M}(t)\right]$ is given
by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)\right]=&C_{d}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\Bigg{[}-\frac{\partial}{\partial
v}\left(I(u,v)\right.\\\
&\left.\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)\right)\Bigg{]}(t-v)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\
+&C_{d}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)du,\\\
\end{aligned}\end{array}$

where $\bar{F}_{\sigma_{M_{s}}}(u)$ and $f_{\sigma_{M_{s}}}(u)$ denote the
survival function and density function of $\sigma_{M_{s}}$ given by (4),
$\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}$ and $I(x,y)$ the survival functions
given by (5) and (6), and $f_{1}(x)$ is the density function of the survival
function given by (7).

For $t\geq T$, the expected transient cost fulfils the following recursive
equation

$E\left[C_{T}^{M}(t)\right]=\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(t-kT)\right]P_{R_{1}}^{M}(kT)+G_{T}^{M}(t),$ (10)

where

$\begin{array}[]{l}\begin{aligned}
G_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\Big{(}C_{p}+C_{I}(k-1)\Big{)}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\Big{(}C_{c}+C_{I}(k-1)\Big{)}P_{R_{1,c}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&\lfloor t/T\rfloor C_{I}\Bigg{(}1-\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)}\\\ +&C_{d}E\left[W_{T}^{M}(\lfloor
t/T\rfloor T,t)\right]\Bigg{(}1-\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)},\end{aligned}\end{array}$

with initial condition $E\left[C_{T}^{M}(0)\right]=0$ and where $\lfloor
t/T\rfloor$ denotes the integer part of $t/T$.

###### Proof.

It is provided in Appendix A. ∎

######  Corollary 1.

Setting $E\left[C_{T}^{M(i)}(t)\right]=E\left[C_{T}^{M}(t)\right]$, for all
$(i-1)T<t\leq iT$ with $i=1,2,\ldots,\lfloor t_{f}/T\rfloor$ the expected
transient cost, $E\left[C_{T}^{M(i)}(t)\right]$, is given by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M(1)}(t)\right]=&C_{d}\Bigg{(}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\Bigg{[}-\frac{\partial}{\partial
v}\left(I(u,v)\right.\\\
&\left.\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)\right)\Bigg{]}(t-v)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\Bigg{)}\\\
+&C_{d}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)du,\\\
\end{aligned}\end{array}$

and for $i\geq 1$

$\begin{array}[]{l}\begin{aligned}
E&\left[C_{T}^{M(i+1)}(t)\right]=G_{T}^{M(i)}(t)\\\
&+\displaystyle\sum_{k=1}^{i}E\left[C_{T}^{M(i+1-k)}(t-kT)\right]P_{R_{1}}^{M}(kT),\end{aligned}\end{array}$

where

$\begin{array}[]{l}\begin{aligned}
G_{T}^{M(i)}(t)=&\displaystyle\sum_{k=1}^{i}\Big{(}C_{p}+C_{I}(k-1)\Big{)}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{i}\Big{(}C_{c}+C_{I}(k-1)\Big{)}P_{R_{1,c}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{i}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&iC_{I}\Big{(}1-\sum_{k=1}^{i}P_{R_{1}}^{M}(kT)\Big{)}\\\
+&C_{d}E\left[W_{T}^{M}(iT,t)\right]\Big{(}1-\sum_{k=1}^{i}P_{R_{1}}^{M}(kT)\Big{)}.\end{aligned}\end{array}$

In order to analyse the uncertainty associated with the expected transient
cost, the standard deviation is calculated. Let
$\Big{(}S^{M}_{T}(t)\Big{)}^{2}$ be the variance of expected transient cost at
time $t$ with periodic inspection times $T$ and preventive threshold $M$
defined as

$\Big{(}S^{M}_{T}(t)\Big{)}^{2}=E\left[C_{T}^{M}(t)^{2}\right]-\left(E\left[C_{T}^{M}(t)\right]\right)^{2}.$
(11)

Based on Theorem 1, the following result is obtained.

######  Theorem 2.

For $t<T$, the expected square cost at time $t>0$,
$E\left[C_{T}^{M}(t)^{2}\right]$, is given by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)^{2}\right]=&C_{d}^{2}\Bigg{(}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\Bigg{[}-\frac{\partial}{\partial
v}\left(I(u,v)\right.\\\
&\left.\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)\right)\Bigg{]}(t-v)^{2}\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\Bigg{)}\\\
+&C_{d}^{2}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)^{2}du.\\\
\end{aligned}\end{array}$

For $t\geq T$, the mean square fulfils the following recursive equation

$E\left[C_{T}^{M}(t)^{2}\right]=\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(t-kT)^{2}\right]P_{R_{1}}^{M}(kT)+H_{T}^{M}(t),$
(12)

where

$\begin{array}[]{l}\begin{aligned}
H_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)^{2}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)+C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\right)^{2}\\\
&P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\\\
&E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,p}}^{M}(kT)\\\ +&\Big{(}\lfloor t/T\rfloor
C_{I}+C_{d}E\left[W_{T}^{M}(\lfloor t/T\rfloor T,t)\right]\Big{)}^{2}\\\
&\Bigg{(}1-\sum_{k=1}^{\lfloor t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)},\\\
\end{aligned}\end{array}$

with initial condition $E\left[C_{T}^{M}(0)^{2}\right]=0$.

###### Proof.

It is given in Appendix B. ∎

######  Corollary 2.

Setting $E\left[C_{T}^{M(i)}(t)^{2}\right]=E\left[C_{T}^{M}(t)^{2}\right]$,
for all $(i-1)T<t\leq iT$ with $i=1,2,\ldots,\lfloor t_{f}/T\rfloor$ the
expected square cost, $E\left[C_{T}^{M(i)}(t)^{2}\right]$, is given by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M(1)}(t)^{2}\right]=&C_{d}^{2}\Bigg{(}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\Bigg{[}-\frac{\partial}{\partial
v}\left(I(u,v)\right.\\\
&\left.\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)\right)\Bigg{]}(t-v)^{2}\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\Bigg{)}\\\
+&C_{d}^{2}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)^{2}du,\\\
\end{aligned}\end{array}$

and for $i\geq 1$

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M(i+1)}(t)^{2}\right]=&\displaystyle\sum_{k=1}^{i}E\left[C_{T}^{M(i+1-k)}(t-kT)^{2}\right]P_{R_{1}}^{M}(kT)\\\
+&H_{T}^{M(i)}(t),\end{aligned}\end{array}$

where

$\begin{array}[]{l}\begin{aligned}
H_{T}^{M(i)}(t)=&\displaystyle\sum_{k=1}^{i}\left(C_{p}+C_{I}(k-1)\right)^{2}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{i}\left(C_{c}+C_{I}(k-1)+C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\right)^{2}\\\
&P_{R_{1,c}}^{M}(kT)\\\
+&2\sum_{k=1}^{i}\left(C_{c}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,c}}^{M}(kT)\\\
+&2\sum_{k=1}^{i}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\\\
&E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&2\sum_{k=1}^{i}\left(C_{p}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,p}}^{M}(kT)\\\
+&\Big{(}iC_{I}+C_{d}E\left[W_{T}^{M}(iT,t)\right]\Big{)}^{2}\\\
&\Bigg{(}1-\sum_{k=1}^{i}P_{R_{1}}^{M}(kT)\Bigg{)}.\\\
\end{aligned}\end{array}$

Hence, by (11) the standard deviation of the transient cost at time $t$,
$S^{M}_{T}(t)$, is given by

$S^{M}_{T}(t)=\sqrt{E\left[C_{T}^{M}(t)^{2}\right]-\left(E\left[C_{T}^{M}(t)\right]\right)^{2}}.$
(13)

##  4 Performance measures of the system

In addition to the expected cost and its standard deviation associated,
recursive expressions for the availability, the reliability and the interval
reliability of the system are obtained.

Let $A_{T}^{M}(t)$ be the availability of the system at time $t>0$, with time
between inspections $T$ and preventive threshold $M$, that is, the probability
that the system is working at time $t$

$\displaystyle A_{T}^{M}(t)$ $\displaystyle=$
$\displaystyle\displaystyle\sum_{j=0}^{\infty}\mathbf{1}_{\\{R_{j}\leq
t<R_{j+1}\\}}P\left[X(t-R_{j})<L,Y>(t-R_{j})\right].$

Next result provides the Markov renewal equation that fulfils the availability
$A_{T}^{M}(t)$.

######  Theorem 3.

For $t<T$, $A_{T}^{M}(t)$ is given by

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\bar{F}_{\sigma_{M_{s}}}(t)\bar{F}_{1}(t)\\\
+&\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-u)I(u,t)du.\\\
\end{aligned}\end{array}$

For $t\geq T$, $A_{T}^{M}(t)$ fulfils the following Markov renewal equation

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}A_{T}^{M}(t-kT)P_{R_{1}}^{M}(kT)\\\
+&J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$
(14)

where

$\begin{array}[]{l}\begin{aligned}
J_{T,1}^{M}(t)=&\bar{F}_{\sigma_{M}}(t)\bar{F}_{1}(t)\\\
+&\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}}(u)\int_{u}^{t}f_{\sigma_{M_{s}}-\sigma_{M}}(v-u)\\\
&\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-v)I(v,t)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\ +&\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}}(u)\bar{F}_{\sigma_{M_{s}}-\sigma_{M}}(t-u)\bar{F}_{1}(t)\leavevmode\nobreak\
du,\\\ \end{aligned}\end{array}$ (15)

and

$\begin{array}[]{l}\begin{aligned}
J_{T,2}^{M}(t)=&\bar{F}_{\sigma_{M_{s}}}(t)\bar{F}_{1}(t)\\\
+&\displaystyle\int_{0}^{\lfloor t/T\rfloor
T}f_{\sigma_{M_{s}}}(u)\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}-\sigma_{M_{s}}}(v-u)\\\
&\bar{F}_{\sigma_{L}-\sigma_{M}}(t-v)I(u,t)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\ +&\displaystyle\int_{0}^{\lfloor t/T\rfloor
T}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{M}-\sigma_{M_{s}}}(t-u)I(u,t)\leavevmode\nobreak\
du\\\ +&\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-u)I(u,t)\leavevmode\nobreak\
du,\\\ \end{aligned}\end{array}$ (16)

with initial condition $A_{T}^{M}(0)=1$ and where $P_{R_{1}}^{M}(kT)$ is given
by (9).

###### Proof.

It is given in Appendix C. ∎

Often, it is also of interest the probability that the system starts working
at time 0 and it continues operating for a time interval. Let $R_{T}^{M}(t)$
be the reliability of the system at time $t$ with time between inspections $T$
and preventive threshold $M$, that is, the probability that the system is
working in $(0,t]$ given by

$\begin{array}[]{l}\begin{aligned} R_{T}^{M}(t)=&P\left[O(u)<L,\forall
u\in(0,t],N_{s}(0,t)=0\right],\end{aligned}\end{array}$

where $O(t)$ denotes the deterioration level of the maintained system at time
$t$, that is,

$O(t)=\sum_{j=0}^{\infty}\mathbf{1}_{\\{R_{j}\leq t<R_{j+1}\\}}X(t-R_{j})$

Based on Theorem 3, the following result is obtained.

######  Theorem 4.

For $t<T$, $R_{T}^{M}(t)$ is given by

$\begin{array}[]{l}\begin{aligned} R_{T}^{M}(t)=&A_{T}^{M}(t).\\\
\end{aligned}\end{array}$

For $t\geq T$, $R_{T}^{M}(t)$ fulfils the following Markov renewal equation

$\begin{array}[]{l}\begin{aligned}
R_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}R_{T}^{M}(t-kT)P_{R_{1,p}}^{M}(kT)\\\
+&J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$
(17)

with initial condition $R_{T}^{M}(0)=1$ where $J_{T,1}^{M}$ and $J_{T,2}^{M}$
are given by (15) and (16), respectively.

###### Proof.

It is given in Appendix D. ∎

A performance measure that extends the reliability is the interval
reliability, defined as the probability that the system is working at time
$t$, and will continue working over a finite time interval of length $s$. The
joint interval reliability is applied when there are periods in the lifetime
cycle when a failure should be avoided with high probability. Let
$IR_{T}^{M}(t,t+s)$ be the interval reliability in $(t,t+s]$. That is

$\displaystyle IR_{T}^{M}(t,t+s)$ $\displaystyle=P\left[O(u)<L,\forall
u\in(t,t+s],N_{s}(t,t+s)=0\right].$

Availability and reliability are particular cases of the interval reliability
since

$A_{T}^{M}(t)=IR_{T}^{M}(t,t+0),$

and

$R_{T}^{M}(t)=IR_{T}^{M}(0,0+t).$

Based on Theorems 3 and 4, the following result is obtained.

######  Theorem 5.

For $t+s<T$, $IR_{T}^{M}(t,t+s)$ is given by

$\begin{array}[]{l}\begin{aligned} IR_{T}^{M}(t,t+s)=&R_{T}^{M}(t+s).\\\
\end{aligned}\end{array}$

For $t+s\geq T$, $IR_{T}^{M}(t,t+s)$ fulfils the following Markov renewal
equation

$\begin{array}[]{l}\begin{aligned}
IR_{T}^{M}(t,t+s)=&\displaystyle\sum_{k=\lfloor
t/T\rfloor+1}^{\lfloor(t+s)/T\rfloor}R_{T}^{M}(t+s-kT)P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}IR_{T}^{M}(t-kT,t+s-kT)P_{R_{1}}^{M}(kT)\\\
+&J_{T,1}^{M}(t+s)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t+s)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$
(18)

with initial conditions $R_{T}^{M}(0)=1$ and $IR_{T}^{M}(0,0)=1$.

###### Proof.

It is given in Appendix E. ∎

##  5 Numerical examples

In this section, some numerical examples are provided to illustrate the
analytical results. To this end, we consider a system subject to an underlying
degradation process modelled as a homogeneous gamma process with parameters
$\alpha=\beta=0.1$. We assume that the system fails when the deterioration
level of the system reaches the breakdown threshold $L=30$. The system can
also fail due to a sudden shock and the sudden shock process is modelled under
a DSPP with intensity

$\lambda\left(t,X(t)\right)=0.01\cdot\mathbf{1}_{\left\\{X(t)\leq
M_{s}\right\\}}+0.1\cdot\mathbf{1}_{\left\\{X(t)>M_{s}\right\\}},\quad t\geq
0,$

where $M_{s}=20$. Under these conditions, the expected time to the failure due
to deterioration is $E\left[\sigma_{L}\right]=34.0335$ $t.u.$ and the expected
time to the failure due to a sudden shock is $E\left[Y\right]=28.3556$ $t.u.$
In addition, we assume the cost sequence $C_{c}=300$ $m.u.$, $C_{p}=150$
$m.u.$, $C_{I}=45$ $m.u.$, and $C_{d}=25$ $m.u./t.u.$ We assume that the life
cycle of the system is $(0,50]$.

The following examples were executed by using MATLAB software version R2014a
on an Intel Core i5-2500 processor with 8GB DDR3 RAM, running under Windows 7
Professional.

###  5.1 Asymptotic expected cost rate analysis

Considering the previous dataset, the expected cost rate is first computed to
establish the values of $T$ and $M$ which will be used subsequently in the
transient approach analysis. In this way, the results obtained by using
asymptotic and transient approaches shall be compared.

The optimisation problem for the expected cost based on the asymptotic formula
given in (8) is computed as follows:

  1. 1.

A grid of size $10$ is obtained discretising the set $[5,50]$ into $10$
equally spaced points from $5$ to $50$ for $T$. For $i=1,2,\ldots,10$, let
$T_{i}$ be the $i$-th value of the grid obtained previously.

  2. 2.

A grid of size $30$ is obtained by discretising the set $[1,30]$ into $30$
equally spaced points from $1$ to $30$ for $M$. For $j=1,2,\ldots,30$, let
$M_{j}$ be the $j$-th value of the grid obtained previously.

  3. 3.

For each fixed combination $(T_{i},M_{j})$, we obtain $50000$ simulations of
$(R_{1},I_{1},W_{d})$, where $R_{1}$ corresponds to the time to the first
replacement (corrective or preventive), $I_{1}$ the nature of the first
maintenance action performed (corrective or preventive), and $W_{d}$ the
downtime up to the first maintenance action.

With these simulations, and applying Monte Carlo method, we obtain
$\tilde{P}_{R_{1,p}}^{M_{j}}(kT_{i})$, $\tilde{P}_{R_{1,c}}^{M_{j}}(kT_{i})$,
$\tilde{P}_{R_{1}}^{M_{j}}(kT_{i})$, and
$\tilde{E}\left[W_{T_{i}}^{M_{j}}((k-1)T_{i},kT_{i}))\right]$ which correspond
to the estimations of $P_{R_{1,p}}^{M_{j}}(kT_{i}))$,
$P_{R_{1,c}}^{M_{j}}(kT_{i}))$, $P_{R_{1}}^{M_{j}}(kT_{i}))$, and
$E\left[W_{T_{i}}^{M_{j}}((k-1)T_{i},kT_{i}))\right]$ for
$k=1,2,\ldots,\lfloor 50/T_{i}\rfloor$, respectively (see Fig. 4).

Fig. 4: Procedure of the Monte Carlo simulation method for variable $T$ and
$M$.

  4. 4.

Quantity $\tilde{C}^{\infty}(T,M)$, which represents the asymptotic expected
cost rate, is calculated by using Equation (8) replacing the corresponding
probabilities by their estimations calculated in Step 3.

  5. 5.

The optimisation problem is reduced to find the values $T_{opt}$ and $M_{opt}$
which minimise the asymptotic expected cost rate $\tilde{C}^{\infty}(T,M)$.
That is

$\tilde{C}^{\infty}(T_{opt},M_{opt})=\min_{\begin{subarray}{c}T\geq 0\\\ 0\leq
M\leq L\end{subarray}}\left\\{\tilde{C}^{\infty}(T,M)\right\\}.$

Fig. 5: Mesh and contour plots for the asymptotic expected cost rate.

Fig. 5 shows the expected cost rate versus $T$ and $M$. The values of $T$ and
$M$ which minimise the expected cost rate are reached at $M_{opt}=14$ $d.u.$
and $T_{opt}=10$ $t.u.$, with an expected cost rate of $15.3819$ $m.u./t.u.$
Below, the expected cost in the finite life cycle will be compared to the
asymptotic expected cost using the values $T_{opt}$ and $M_{opt}$.

###  5.2 Expected transient cost rate analysis for a fixed $T$

We consider a time between inspections $T=10$ $t.u.$ The optimisation problem
for the expected transient cost based on the recursive formula given in (10)
is computed as follows:

  1. 1.

A grid of size $30$ is obtained by discretising the set $[1,30]$ into $30$
equally spaced points from $1$ to $30$ for $M$.

  2. 2.

For fixed $T=10$ $t.u.$ and for each fixed $M_{j}$, we obtain $50000$
simulations of $(R_{1},I_{1},W_{d})$. With these simulations, and applying
Monte Carlo method, we obtain the estimations
$\tilde{P}_{R_{1,p}}^{M_{j}}(10k)$, $\tilde{P}_{R_{1,c}}^{M_{j}}(10k)$,
$\tilde{P}_{R_{1}}^{M_{j}}(10k)$, and
$\tilde{E}\left[W_{10}^{M_{j}}((k-1)10,10k)\right]$ (see Fig. 4).

  3. 3.

For fixed $T=10$ $t.u.$ and $M_{j}$, let
$\tilde{E}\left[C_{10}^{M_{j}}(50)\right]$ be the expected cost in the finite
life cycle. The expected transient cost is calculated by using the recursive
formula given in (10), replacing the corresponding probabilities by their
estimations calculated in Step 2, with initial condition
$\tilde{E}\left[C_{10}^{M_{j}}(0)\right]=0$.

  4. 4.

For fixed $T=10$ $t.u.$, the optimisation problem is reduced to find the value
$M_{opt}$ which minimises the expected cost
$\tilde{E}\left[C_{10}^{M}(50)\right]$. That is

$\tilde{E}\left[C_{10}^{M_{opt}}(50)\right]=\min_{0\leq M\leq
L}\left\\{\tilde{E}\left[C_{10}^{M}(50)\right]\right\\}.$

Let $\tilde{E}\left[C_{10}^{M_{j}}(50)\right]/50$ be the expected cost rate in
the life cycle of the system. The expected cost rate calculated using the
recursive method and the expected cost rate calculated using strictly Monte
Carlo simulation are shown in Fig. 6. The expected transient cost rate based
on strictly Monte Carlo simulation was calculated for $30$ equally spaced
points in the interval $(0,30]$ with $50000$ simulations for each point. Based
on Fig. 6, for the recursive method, the expected transient cost rate based on
the recursive method given by Theorem 1 reaches its minimum value at $M=14$
$d.u.$, with an expected transient cost rate of $14.7639$ $m.u./t.u.$ On the
other hand, the expected transient cost rate based on strictly Monte Carlo
simulation reaches its minimum value at $M=14$ $d.u.$, with an expected
transient cost rate of $15.2096$ $m.u./t.u.$

Fig. 6: Expected cost rate in the life cycle for different values of $M$.

For $T=10$ $t.u.$, Fig. 7 shows the expected transient cost rate calculated by
using the recursive method and the asymptotic expected cost rate calculated
throughout the procedure detailed in Section 5.1 versus $M$.

Fig. 7: Expected transient cost rate and asymptotic expected cost rate at time
$t_{f}=50$ $t.u.$ for different values of $M$.

As we said previously, the value of $M$ which minimises the expected transient
cost rate is reached at $M=14$ $d.u.$, with an expected cost rate of $14.7639$
$m.u./t.u.$ On the other hand, the asymptotic expected cost rate reaches its
minimum value at $M=14$ $d.u.$, with an expected cost rate of $15.3819$
$m.u./t.u.$

Now, we calculate the standard deviation of the cost. Fig. 8 shows the
expected transient cost rate with its standard deviation associated given by
Equation (13). Both quatities were calculated for $30$ equally spaced points
in $(0,30]$ by using the recursive formula given in (12) throughout the steps
detailed in 5.2.

Fig. 8: Expected cost rate in the life cycle and standard deviation versus
$M$.

Now, we focus on the influence of the main model parameters on the expected
cost in the life cycle. Firstly, a sensitivity analysis of the gamma process
parameters is performed.

The values of the gamma process parameters are modified according to the
following specifications:

$\alpha_{(v_{i}\%)}=\alpha\left[1+\frac{v_{i}}{100}\right]\quad\text{and}\quad\beta_{(v_{j}\%)}=\beta\left[1+\frac{v_{j}}{100}\right],$
(19)

where $v_{i}$ and $v_{j}$ are, respectively, the $i$-th and $j$-th position of
the vector $\mathbf{v}=(-10$, $-5$, $-1$, $0$, $1$, $5$, $10$). Then, the
parameter values for $\alpha$ and $\beta$ can be simultaneous and
independently modified both for increasing and decreasing changes.

Let
$\tilde{E}\left[C_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(t_{f})\right]$
be the minimal expected transient cost at time $t_{f}$ obtained when the gamma
process parameters ($\alpha$ and $\beta$) are varied according to the
specifications given in (19). The expected transient cost for each combination
of $\alpha_{(v_{i}\%)}$ and $\beta_{(v_{j}\%)}$ are calculated based on the
recursive method following the steps detailed in 5.2. The relative measure
$V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$ is defined as

$\frac{\left|\tilde{E}\left[C_{10}^{M_{opt}}(50)\right]-\tilde{E}\left[C_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)\right]\right|}{\tilde{E}\left[C_{10}^{M_{opt}}(50)\right]},$
(20)

where $\tilde{E}\left[C_{10}^{M_{opt}}(50)\right]$ is the minimal expected
transient cost calculated in 5.2.

For fixed $i$ and $j$, $V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$
measures the relative difference between the minimal expected transient cost
with the original parameter values and the minimal expected transient cost
calculated using the modified parameter values for fixed $M$ and $T=10$ $t.u.$
Values closer to zero have a lower influence on the expected transient cost
rate.

Table 4 shows the relative variation percentages with a shaded grey scale.
Each cell represents $V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$
expressed in percentage. Darker colours of cells denote a higher relative
variation percentage. The results obtained show that
$V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$ grows when $\alpha$
increases and $\beta$ decreases and
$V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$ decreases when $\alpha$
decreases and $\beta$ increases. In this way,
$V_{10,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{M}(50)$ reaches its minimum
value when $\alpha$ is minimum and $\beta$ is maximum and its maximum value
when $\alpha$ is maximum and $\beta$ is minimum.

| $\beta_{{}_{(-10\%)}}$ | $\beta_{{}_{(-5\%)}}$ | $\beta_{{}_{(-1\%)}}$ | $\beta$ | $\beta_{{}_{(1\%)}}$ | $\beta_{{}_{(5\%)}}$ | $\beta_{{}_{(10\%)}}$
---|---|---|---|---|---|---|---
$\alpha_{(-10\%)}$ | 1.1862 | 2.2865 | 5.8303 | 5.7383 | 6.4372 | 8.6707 | 11.0707
$\alpha_{(-5\%)}$ | 4.1067 | 0.9778 | 1.9635 | 2.7403 | 3.2077 | 6.0850 | 9.0580
$\alpha_{(-1\%)}$ | 7.0083 | 2.7336 | 0.6373 | 0.3666 | 1.0904 | 3.6113 | 6.1380
$\alpha$ | 8.0187 | 3.8252 | 1.0523 | 0.0000 | 0.2624 | 2.8032 | 6.0328
$\alpha_{(1\%)}$ | 8.7105 | 4.5726 | 1.4507 | 1.0696 | 0.0267 | 2.6336 | 5.1357
$\alpha_{(5\%)}$ | 11.3000 | 7.0853 | 4.0552 | 3.6551 | 2.7247 | 0.2344 | 2.7598
$\alpha_{(10\%)}$ | 15.2929 | 10.5769 | 7.2532 | 6.6684 | 5.7666 | 3.0604 | 0.0958

Table 1: Relative variation percentages for the expected transient cost for
the gamma process parameters for a fixed $T=10$ $t.u.$

By modifying $\pm 1\%$ around $\alpha=\beta=0.1$, the relative variation
percentages are small. The results also show that the relative variation
percentages are lower in the diagonal of the table. That means when the
parameters $\alpha$ and $\beta$ are modified in the same direction and
magnitude.

Similarly, the values of the parameters $\lambda_{1}$ and $\lambda_{2}$ are
modified according to the following specifications:

$\lambda_{1,(v_{i}\%)}=\lambda_{1}\left[1+\frac{v_{i}}{100}\right]\quad\text{and}\quad\lambda_{2,(v_{j}\%)}=\lambda_{2}\left[1+\frac{v_{j}}{100}\right],$
(21)

Let
$\tilde{E}^{*}\left[C^{M}_{10,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(t_{f})\right]$
be the minimal expected transient cost obtained by varying the parameters
$\lambda_{1}$ and $\lambda_{2}$ simultaneously as in the scheme given in (21).
Now, the relative variation
$V^{M}_{10,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(50)$ is given by

$\frac{\left|\tilde{E}^{*}\left[C_{10}^{M_{opt}}(50)\right]-\tilde{E}^{*}\left[C_{10,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}^{M}(50)\right]\right|}{\tilde{E}^{*}\left[C_{10}^{M_{opt}}(50)\right]}.$
(22)

The relative variation percentages are presented in Table 5. The results show
that the parameter $\lambda_{1}$ has greater effects on
$V^{M}_{10,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(50)$ than the
parameter $\lambda_{2}$, reaching the lowest values when the variation for
$\lambda_{1}=0.01$ is minimal, that is $\pm 1\%$, and the highest values when
the variation for $\lambda_{1}$ is maximised, that is $\pm 10\%$.

| $\lambda_{2,(-10\%)}$ | $\lambda_{2,(-5\%)}$ | $\lambda_{2,(-1\%)}$ | $\lambda_{2}$ | $\lambda_{2,(1\%)}$ | $\lambda_{2,(5\%)}$ | $\lambda_{2,(10\%)}$
---|---|---|---|---|---|---|---
$\lambda_{1,(-10\%)}$ | 3.0794 | 2.4761 | 2.3266 | 1.9374 | 1.8849 | 2.0284 | 1.3703
$\lambda_{1,(-5\%)}$ | 1.6095 | 1.2419 | 1.3581 | 0.8946 | 0.9002 | 0.3058 | 0.3755
$\lambda_{1,(-1\%)}$ | 0.6159 | 0.5227 | 0.2235 | 0.0000 | 0.5353 | 0.0241 | 0.9365
$\lambda_{1}$ | 0.4598 | 0.2900 | 0.2235 | 0.0000 | 0.5353 | 0.0241 | 0.9365
$\lambda_{1,(1\%)}$ | 0.4598 | 0.2451 | 0.4946 | 0.0118 | 0.3289 | 1.0549 | 0.9833
$\lambda_{1,(5\%)}$ | 0.2698 | 0.3685 | 0.5846 | 0.7321 | 0.4718 | 0.8365 | 1.1554
$\lambda_{1,(10\%)}$ | 1.3381 | 2.3405 | 2.7045 | 2.2538 | 2.7964 | 2.5593 | 3.2779

Table 2: Relative variation percentages for the expected transient cost for
parameters $\lambda_{1}$ and $\lambda_{2}$ for a fixed $T=10$ $t.u.$

###  5.3 Analysis of the expected cost rate in the life cycle for a fixed $M$

We now analyse the time between inspections $T$ influence on the expected cost
in the life cycle of the system for a preventive threshold $M=14$ $d.u.$ As in
Section 5.2, the optimisation problem for the expected transient cost based on
the recursive formula given in (10) is computed throughout the following
steps:

  1. 1.

A grid of size $10$ is obtained by discretising the set $[5,50]$ into $10$
equally spaced points from $5$ to $50$ for the time between inspections $T$.

  2. 2.

For fixed $M=14$ $d.u.$ and for each fixed $T_{i}$, we obtain $50000$
simulations of $(R_{1},I_{1},W_{d})$. With these simulations and applying
Monte Carlo method, we obtain the estimations
$\tilde{P}_{R_{1,p}}^{14}(kT_{i})$, $\tilde{P}_{R_{1,c}}^{14}(kT_{i})$,
$\tilde{P}_{R_{1}}^{14}(kT_{i})$, and
$\tilde{E}\left[W_{T_{i}}^{14}((k-1)T_{i},kT_{i})\right]$ for
$k=1,2,\ldots,\lfloor 50/T_{i}\rfloor$ (see Fig. 4).

  3. 3.

Let $\tilde{E}\left[C_{T_{i}}^{14}(50)\right]$ be the expected cost in the
life cycle. The expected cost is calculated by using the recursive formula
given in (10), replacing the corresponding probabilities by their estimations
calculated in Step 2, with initial condition
$\tilde{E}\left[C_{T_{i}}^{14}(0)\right]=0$.

  4. 4.

For fixed $M=14$ $d.u.$, the optimisation problem is reduced to find the value
of $T_{opt}$ which minimises $\tilde{E}\left[C_{T}^{14}(50)\right]$, that is,

$\tilde{E}\left[C_{T_{opt}}^{14}(50)\right]=\min_{T\geq
0}\left\\{\tilde{E}\left[C_{T}^{14}(50)\right]\right\\}.$

Let $\tilde{E}\left[C_{T_{i}}^{14}(50)\right]/50$ be the expected cost rate in
the life cycle of the system. The expected cost rate evaluated using the
recursive formula and using Monte Carlo simulation are shown in Fig. 9. The
expected cost rate based on strictly Monte Carlo simulation was calculated for
$10$ equally spaced points in the interval $(0,50]$ with $50000$ realizations
for each point.

Fig. 9: Expected transient cost rate at time $t_{f}=50$ $t.u.$ versus $T$.

Based on Fig. 9, for the recursive method, the expected cost rate at time
$t_{f}=50$ $t.u.$ reaches its minimum value for $T=10$ $t.u.$, with an
expected cost rate of $14.7637$ $m.u./t.u.$ On the other hand, using strictly
Monte Carlo simulation, the expected cost rate at time $t_{f}=50$ $t.u.$
reaches its minimum value at $T=30$ $t.u.$, with an expected cost rate of
$12.8252$ $m.u./t.u.$ This difference can be explained due to the high
variance in deterioration increments of the degradation process.

For $M=14$ $d.u.$ Table 3 shows the average number of renewals of the system
in its life cycle for each value of $T$.

$\mathbf{T}$ | 5 | 10 | 15 | 20 | 25
---|---|---|---|---|---
$\mathbf{E\left[N_{T}^{14}(50)\right]}$ | 2.4650 | 2.2007 | 1.7772 | 1.4327 | 1.6419
$\mathbf{T}$ | 30 | 35 | 40 | 45 | 50
$\mathbf{E\left[N_{T}^{14}(50)\right]}$ | 0.8884 | 0.93964 | 0.9682 | 0.9833 | 0.9926

Table 3: Average number of complete renewal cycles up to $t_{f}=50$ for
different values of $T$.

For fixed $M=14$ $d.u.$, Fig. 10 shows the expected cost rate in the life
cycle of the system calculated by using the recursive method and the
asymptotic expected cost rate calculated throughout the procedure detailed in
Section 5.1 versus $T$.

Fig. 10: Expected transient cost rate and asymptotic expected cost rate at
time $t_{f}=50$ $t.u.$ for different values of $T$.

As we said previously, the expected transient cost rate calculated using the
recursive method at time $t_{f}=50$ $t.u.$ reaches its minimum value at $T=10$
$t.u.$, with an expected transient cost rate of $14.7637$ $m.u./t.u.$ On the
other hand, the asymptotic expected cost rate reaches its minimum value at
$T=10$ $t.u.$, with an asymptotic expected cost rate of $15.3819$ $m.u./t.u.$
Asymptotic expected cost rate shows a smoother behaviour compared to the
expected transient cost rate.

Next, the standard deviation of the cost is calculated. Fig. 11 shows the
expected transient cost rate calculated using the recursive method with its
standard deviation associated given in (13) versus $T$. This deviation was
calculated for $10$ equally spaced points in the interval $(0,50]$ by using
the recursive formula given in Equation (12) throughout the steps detailed in
5.5.

Fig. 11: Expected transient cost rate at time $t_{f}=50$ $t.u.$ and its
standard deviation associated for different values of $T$.

Focusing now on the main model parameters influence on the solution, we
analyse the gamma process parameters sensitivity.

Let
$\tilde{E}\left[C_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{14}(50)\right]$ be
the minimal expected transient cost at time $t_{f}=50$ $t.u.$ obtained when
the gamma process parameters ($\alpha$ and $\beta$) are varied according to
the specifications given in (19) for each value of $T$ with preventive
threshold $M=14$ $d.u.$ Based on (20),
$V_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{14}(50)$ denotes the relative
variation between the minimal expected transient cost with the original
parameter values and the minimal expected transient cost calculated by using
the parameter values modified according to (19) for variable $T$ and $M=14$
$d.u.$

Table 4 shows the values obtained for
$V_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{14}(50)$ expressed in percentage.
By modifying $\pm 1\%$ around $\alpha=\beta=0.1$, the relative variation
percentages are small. The results obtained also show that
$V_{T,\alpha_{(v_{i}\%)},\beta_{(v_{j}\%)}}^{14}(50)$ is lower in the diagonal
of the table. That means when the parameters $\alpha$ and $\beta$ are modified
in the same direction and magnitude.

| $\beta_{{}_{(-10\%)}}$ | $\beta_{{}_{(-5\%)}}$ | $\beta_{{}_{(-1\%)}}$ | $\beta$ | $\beta_{{}_{(1\%)}}$ | $\beta_{{}_{(5\%)}}$ | $\beta_{{}_{(10\%)}}$
---|---|---|---|---|---|---|---
$\alpha_{(-10\%)}$ | 1.5684 | 1.8862 | 4.7069 | 5.7107 | 6.1450 | 8.0609 | 11.2240
$\alpha_{(-5\%)}$ | 4.5915 | 0.9400 | 1.7235 | 2.5432 | 3.3861 | 5.1360 | 8.1300
$\alpha_{(-1\%)}$ | 8.5189 | 3.4618 | 0.9003 | 0.0450 | 0.6994 | 3.0001 | 6.2257
$\alpha_{1}$ | 4.4937 | 8.2817 | 1.3744 | 0.0000 | 0.1749 | 2.7793 | 5.2119
$\alpha_{(1\%)}$ | 8.9568 | 4.7309 | 1.9896 | 1.8475 | 0.4936 | 1.4929 | 5.0625
$\alpha_{(5\%)}$ | 11.5044 | 7.2387 | 5.0048 | 3.5146 | 3.3441 | 0.8091 | 2.7099
$\alpha_{(10\%)}$ | 15.1620 | 10.7048 | 7.5764 | 6.5529 | 6.7841 | 3.6827 | 0.6109

Table 4: Relative variation percentages for the expected transient cost for
the gamma process parameters for a fixed $M=14$ $d.u.$

On the other hand, let
$\tilde{E}^{*}\left[C^{14}_{T,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(t_{f})\right]$
be the minimal expected transient cost obtained by varying the parameters
$\lambda_{1}$ and $\lambda_{2}$ simultaneously as in the scheme given in (21).
Based on (22), $V^{14}_{T,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(50)$
denotes the relative variation between the minimal expected transient cost
with the original parameter values and the minimal expected transient cost
calculated by using the parameter values modified according to (21) for
variable $T$ and $M=14$ $d.u.$

The relative variation percentages are presented in Table 5. The results show
that when $\lambda_{1}=0.01$ is modified between $-5\%$ and $1\%$, the
relative variation percentages are small. In addition, the parameter
$\lambda_{1}$ has greater effects on
$V^{14}_{T,\lambda_{1,(v_{i}\%)},\lambda_{2,(v_{j}\%)}}(50)$ than the
parameter $\lambda_{2}$.

| $\lambda_{2,(-10\%)}$ | $\lambda_{2,(-5\%)}$ | $\lambda_{2,(-1\%)}$ | $\lambda_{2}$ | $\lambda_{2,(1\%)}$ | $\lambda_{2,(5\%)}$ | $\lambda_{2,(10\%)}$
---|---|---|---|---|---|---|---
$\lambda_{1,(-10\%)}$ | 2.1506 | 1.9893 | 1.9244 | 1.2564 | 1.7375 | 1.4737 | 1.7043
$\lambda_{1,(-5\%)}$ | 0.8525 | 0.3948 | 0.7135 | 0.2861 | 0.0258 | 0.2199 | 0.2733
$\lambda_{1,(-1\%)}$ | 0.8510 | 0.5039 | 0.4442 | 0.6259 | 0.4789 | 0.7896 | 1.3432
$\lambda_{1}$ | 0.4986 | 0.5000 | 0.4663 | 0.0000 | 0.7516 | 1.0475 | 1.1239
$\lambda_{1,(1\%)}$ | 0.0194 | 0.6036 | 0.9957 | 0.9308 | 0.7760 | 1.1128 | 1.9800
$\lambda_{1,(5\%)}$ | 1.9045 | 1.5869 | 2.8829 | 1.4909 | 2.1346 | 2.5697 | 3.0034
$\lambda_{1,(10\%)}$ | 2.4267 | 2.1395 | 2.8044 | 2.9864 | 2.8605 | 3.3628 | 3.0087

Table 5: Relative variation percentages for the expected transient cost rate
for parameters $\lambda_{1}$ and $\lambda_{2}$ for a fixed $M=14$ $d.u.$

###  5.4 Transient two-dimensional expected cost rate analysis

The expected transient cost based on the recursive formula given in (10)
considering $M$ and $T$ variables is analysed. The optimisation problem is
computed as follows:

  1. 1.

A grid of size $10$ is obtained by discretising the set $[5,50]$ into $10$
equally spaced points from $5$ to $50$ for $T$. Let $T_{i}$ be the $i$-th
value of the grid obtained previously, for $i=1,2,\ldots,10$.

  2. 2.

A grid of size $30$ is obtained by discretising the set $[1,30]$ into $30$
equally spaced points from $1$ to $30$ for $M$. Let $M_{j}$ be the $j$-th
value of $M$ which corresponds to the $i$-th value of the grid obtained
previously, for $j=1,2,\ldots,30$.

  3. 3.

For each fixed combination $(T_{i},M_{j})$, we obtain $50000$ simulations of
$(R_{1},I_{1},W_{d})$. With these simulations and applying Monte Carlo method,
$\tilde{P}_{R_{1,p}}^{M_{j}}(kT_{i})$, $\tilde{P}_{R_{1,c}}^{M_{j}}(kT_{i})$,
$\tilde{P}_{R_{1}}^{M_{j}}(kT_{i})$, and
$\tilde{E}\left[W_{T_{i}}^{M_{j}}((k-1)T_{i},kT_{i}))\right]$ for
$k=1,2,\ldots,\lfloor 50/T_{i}\rfloor$ are obtained (see Fig. 4).

  4. 4.

For fixed $T_{i}$ and $M_{j}$, let
$\tilde{E}\left[C_{T_{i}}^{M_{j}}(50)\right]$ be the expected cost at time
$t_{f}=50$ $t.u.$ Expected cost is calculated by using the recursive formula
given in (10), replacing the corresponding probabilities by their estimations
calculated in Step 2, with initial condition
$\tilde{E}\left[C_{T_{i}}^{M_{j}}(0)\right]=0$.

  5. 5.

The optimisation problem is reduced to find the values $T_{opt}$ and $M_{opt}$
which minimise the expected cost
$\tilde{E}\left[C_{T_{i}}^{M_{j}}(50)\right]$. That is

$\tilde{E}\left[C_{T_{opt}}^{M_{opt}}(50)\right]=\min_{\begin{subarray}{c}T\geq
0\\\ 0\leq M\leq
L\end{subarray}}\left\\{\tilde{E}\left[C_{T}^{M}(50)\right]\right\\}.$

Let $\tilde{E}\left[C_{T_{i}}^{M_{j}}(50)\right]/50$ be the expected cost rate
at time $t_{f}=50$ $t.u.$ The expected cost rate versus $T$ and $M$ is shown
in Fig. 12.

Fig. 12: Mesh and contour plots for the expected transient cost rate.

The values of $T$ and $M$ which minimise the expected cost rate at time
$t_{f}=50$ $t.u.$ are reached for $M=14$ $d.u.$ and $T=10$ $t.u.$ with an
expected cost rate of $14.7637$ $m.u./t.u.$

###  5.5 Performance measures for optimal values of $T$ and $M$

We now analyse the availability and reliability of the system considering the
optimal maintenance strategy.

The availability of the system based on the recursive formula given in (14) is
computed throughout the following steps:

  1. 1.

A grid of size $50$ is obtained by discretising the set $[1,50]$ into $50$
equally spaced points from $1$ to $50$ for the instant time $t$. Let $t_{n}$
be the $n$-th value of the grid obtained previously, for $n=1,2,\ldots,50$.

  2. 2.

Let $\tilde{A}_{10}^{14}(t_{n})$ be the system availability estimation for a
time between inspections $T_{opt}=10$ $t.u.$ and a preventive threshold
$M_{opt}=14$ $d.u.$ The availability of the system is calculated by using the
recursive formula given in (14), replacing $P_{R_{1}}^{14}(10k)$ by its
estimation $\tilde{P}_{R_{1}}^{14}(10k)$ calculated in Step 2 of the procedure
detailed in Section 5.4, with initial condition $\tilde{A}_{10}^{14}(0)=1$.

Figure 13 shows the availability of the system versus $t$. We can conclude
that, for fixed $T=10$ and $M=14$, the probability that the system is working
at any instant time of its life cycle is, at least, of the $82\%$.

Fig. 13: Availability of the system for different values of $t$.

Next, the reliability of the system is evaluated. The reliability of the
system based on the recursive formula given in (17) is computed throughout the
following steps:

  1. 1.

A grid of size $50$ is obtained by discretising the set $(1,50]$ into $50$
equally spaced points from $1$ to $50$ for $t$.

  2. 2.

For fixed $T=10$ and $M=14$, let $\tilde{R}_{10}^{14}(t)$ be the system
reliability estimation at time $t$. The reliability of the system is
calculated by using the recursive formula given in (17), replacing
$P_{R_{1,p}}^{14}(10k)$ by its estimation $\tilde{P}_{R_{1,p}}^{14}(10k)$
calculated in Step 2 of the procedure detailed in Section 5.4 with initial
condition $\tilde{R}_{10}^{14}(0)=1$.

Figure 14 shows the reliability of the system versus $t$. We can conclude
that, for fixed $T=10$ and $M=14$, the probability that the system does not
fail in its life cycle is of the $32\%$.

Fig. 14: Reliability of the system for different values of $t$.

Finally, the interval reliability of the system based on (18) is computed
throughout the following steps:

  1. 1.

A grid of size $10$ is obtained by discretising the set $[10,30]$ into $10$
equally spaced points from $10$ to $30$ for $t$.

  2. 2.

For $s=5$, let $\tilde{IR}_{10}^{14}(t,t+s)$ be the system interval
reliability estimation for a time between inspections $T=10$ $t.u.$ and a
preventive threshold $M=14$ $d.u.$. The interval reliability of the system is
calculated by using the recursive formula given in (18), replacing
$P_{R_{1,p}}^{14}(10k)$ and $P_{R_{1}}^{14}(10k)$ by the estimation
$\tilde{P}_{R_{1,p}}^{14}(10k)$ and $\tilde{P}_{R_{1}}^{14}(10k)$,
respectively, calculated in Step 2 of the procedure detailed in Section 5.4
with initial conditions $\tilde{IR}_{10}^{14}(0,0)=1$ and
$\tilde{R}_{10}^{14}(0)=1$.

Figure 15 shows the interval reliability of the system versus $t$. As we can
observe, the results provide for both methods are very similar. Furthermore we
can conclude that, for fixed $T=10$ and $M=14$, the probability that the
system does not fail in the interval $[t,t+5]$ is, at least, of the $72\%$ for
$15\leq t\leq 35$.

Fig. 15: Interval reliability of the system for different values of $t$.

##  6 Conclusions and further works

In this paper, a CBM strategy is analysed by considering a finite life cycle
of the system. The system is subject to two different causes of failure, a
degradation process modelled under a gamma process and a sudden shock process
which follows a DSPP. We consider that both causes of failure are dependent.
This dependence is reflected in that the system is more susceptible to
external shocks when the deterioration level of the system reaches a certain
threshold $M_{s}$.

Under these assumptions, the expected cost rate in the life cycle is used as
objective function to obtain the optimal maintenance strategy. To this end, a
numerical method based on a recursive formula is provided to evaluate the
expected cost rate and the standard deviation associated.

The expected transient cost rate calculated using the recursive formula is
compared to both the asymptotic expected cost rate and the expected transient
cost rate calculated using strictly Monte Carlo simulation. In addition, the
robustness of the gamma process parameters and DSPP is analysed. Furthermore,
for the comparison between transient and asymptotic cost rate the results are
also provided under a bivariate case, where the time between inspections and
the preventive threshold vary simultaneously. In the comparison between the
method based on strictly Monte Carlo simulation and the method based on the
recursive formula, we observed similar results for a fixed time between
inspections $T$. For the preventive threshold $M$, the results presented some
differences, which could be explained due to the high variance in the
deterioration increments of the degradation process. If the life cycle
increases, the recursive method shall tends to be more costly in terms of
computation than the method based on strictly Monte Carlo simulation.

Finally, three important performance measures in the maintenance field, the
availability, the reliability and the interval reliability of the system, are
analysed under the bivariate optimal maintenance strategy. Recursive formulas
are given to obtain these performance measures. Furthermore, the results
obtained using the recursive formulation is compared to the results obtained
using strictly Monte Carlo simulation.

In this paper, we consider a system subject to a unique degradation process.
However, sometimes the system is subject to multiple degradation processes. A
possible further extension of this work is to consider a system subject to
multiple degradation processes. With respect to the two types of failure, the
analysis of this model is based on the dependence of the degradation level of
the system on the intensity of shocks. An interesting extension could be to
assume a bidirectional relation of dependence where the process affects also
to the degradation process.

## Appendix A

For $t<T$, $E\left[C_{T}^{M}(t)\right]$, is given by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)\right]=&C_{d}E\left[(t-Y)\mathbf{1}_{\\{\sigma_{M_{s}}<Y<t,\leavevmode\nobreak\
Y<\sigma_{L}\\}}\right]\\\
+&C_{d}E\left[(t-\sigma_{L})\mathbf{1}_{\\{\sigma_{M_{s}}<\sigma_{L}<t,\leavevmode\nobreak\
\sigma_{L}<Y\\}}\right]\\\
+&C_{d}E\left[(t-Y)\mathbf{1}_{\\{Y<t,\leavevmode\nobreak\
Y<\sigma_{M_{s}}\\}}\right].\\\ \end{aligned}\end{array}$

That is

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)\right]=&C_{d}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\left[-\frac{\partial}{\partial
v}I(u,v)\right]\\\
&\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)(t-v)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\
+&C_{d}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}I(u,v)\\\
&f_{\sigma_{L}-\sigma_{M_{s}}}(v-u)(t-v)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\
+&C_{d}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)du.\\\
\end{aligned}\end{array}$

For $t\geq T$, $E\left[C_{T}^{M}(t)\right]$ is conditioned to $R_{1}$

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)\right]=&E\left[C_{T}^{M}(t),R_{1}\leq
t\right]+E\left[C_{T}^{M}(t),R_{1}>t\right].\end{aligned}\end{array}$

Thus, if $R_{1}>t$

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(t),R_{1}>t\right]\\\
=&\lfloor t/T\rfloor C_{I}\Bigg{(}1-\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)}\\\ +&C_{d}E\left[W_{T}^{M}(\lfloor
t/T\rfloor T,t)\right]\Bigg{(}1-\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)}.\end{aligned}\end{array}$

If $R_{1}\leq t$, $E\left[C_{T}^{M}(t)\right]$ can be split into two terms:
the cost in the first renewal cycle $\left(C_{T}^{M}(R_{1})\right)$ and the
cost in the remaining time horizon $\left(C_{T}^{M}(R_{1},t)\right)$. Since
$C_{T}^{M}(R_{1})$ and $C_{T}^{M}(R_{1},t)$ are independent, we get

$\begin{array}[]{l}\begin{aligned} E\left[C_{T}^{M}(t),R_{1}\leq
t\right]=&E\left[C_{T}^{M}(R_{1}),R_{1}\leq t\right]\\\
+&E\left[C_{T}^{M}(R_{1},t),R_{1}\leq t\right].\end{aligned}\end{array}$

For a fixed $T$

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(R_{1}),R_{1}\leq
t\right]\\\ =&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(R_{1}),R_{1,c}=kT\right]\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(R_{1}),R_{1,p}=kT\right],\end{aligned}\end{array}$

being

$\begin{array}[]{l}\begin{aligned}
E&\left[C_{T}^{M}(R_{1}),R_{1,c}=kT\right]\\\
=&\left(C_{c}+C_{I}(k-1)\right)P_{R_{1,c}}^{M}(kT)\\\
+&C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]P_{R_{1,c}}^{M}(kT),\end{aligned}\end{array}$

and

$E\left[C_{T}^{M}(R_{1}),R_{1,p}=kT\right]=\left(C_{p}+C_{I}(k-1)\right)P_{R_{1,p}}^{M}(kT).$

Since $C_{T}^{M}(R_{1},t)$ is stochastically the same as $C_{T}^{M}(t-R_{1})$,

$E\left[C_{T}^{M}(R_{1},t),R_{1}=kT\right]=E\left[C_{T}^{M}(t-kT)\right]P_{R_{1}}^{M}(kT).$

Hence, $E\left[C_{T}^{M}(t)\right]$ verifies the following recursive equation

$E\left[C_{T}^{M}(t)\right]=\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(t-kT)\right]P_{R_{1}}^{M}(kT)+G_{T}^{M}(t),$

being

$\begin{array}[]{l}\begin{aligned}
G_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\Big{(}C_{p}+C_{I}(k-1)\Big{)}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\Big{(}C_{c}+C_{I}(k-1)\Big{)}P_{R_{1,c}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&\lfloor t/T\rfloor C_{I}\Big{(}1-\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Big{)}\\\ +&C_{d}E\left[W_{T}^{M}(\lfloor
t/T\rfloor T,t)\right]\Big{(}1-\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Big{)},\\\ \end{aligned}\end{array}$

and the result holds.

## Appendix B

For $t<T$, the expected square cost, $E\left[C_{T}^{M}(t)^{2}\right]$, is
given by

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)^{2}\right]=&C_{d}^{2}E\left[(t-Y)^{2}\mathbf{1}_{\\{\sigma_{M_{s}}<Y<t,\leavevmode\nobreak\
Y<\sigma_{L}\\}}\right]\\\
+&C_{d}^{2}E\left[(t-\sigma_{L})^{2}\mathbf{1}_{\\{\sigma_{M_{s}}<\sigma_{L}<t,\leavevmode\nobreak\
\sigma_{L}<Y\\}}\right]\\\
+&C_{d}^{2}E\left[(t-Y)^{2}\mathbf{1}_{\\{Y<t,\leavevmode\nobreak\
Y<\sigma_{M_{s}}\\}}\right].\\\ \end{aligned}\end{array}$

That is

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)\right]=&C_{d}^{2}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}\left[-\frac{\partial}{\partial
v}I(u,v)\right]\\\
&\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(v-u)(t-v)^{2}\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\
+&C_{d}^{2}\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\int_{u}^{t}I(u,v)\\\
&f_{\sigma_{L}-\sigma_{M_{s}}}(v-u)(t-v)^{2}\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\
+&C_{d}^{2}\displaystyle\int_{0}^{t}f_{1}(u)\bar{F}_{\sigma_{M_{s}}}(u)(t-u)^{2}du.\\\
\end{aligned}\end{array}$

For $t\geq T$, $E\Big{[}C_{T}^{M}(t)^{2}\Big{]}$ is conditioned to $R_{1}$

$\begin{array}[]{l}\begin{aligned}
E\left[C_{T}^{M}(t)^{2}\right]=&E\left[C_{T}^{M}(t)^{2},R_{1}\leq
t\right]+E\left[C_{T}^{M}(t)^{2},R_{1}>t\right].\end{aligned}\end{array}$

Hence,

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(t)^{2},R_{1}>t\right]\\\
=&\Big{(}\lfloor t/T\rfloor C_{I}+C_{d}E\left[W_{T}^{M}(\lfloor t/T\rfloor
T,t)\right]\Big{)}^{2}\\\ &\Bigg{(}1-\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)},\end{aligned}\end{array}$

where $P_{R_{1}}^{M}$ is given by (9). On the other hand

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(t)^{2},R_{1}\leq
t\right]\\\
=&E\left[\left(C_{T}^{M}(R_{1})+C_{T}^{M}(R_{1},t)\right)^{2},R_{1}\leq
t\right].\end{aligned}\end{array}$

Since $C_{T}^{M}(R_{1})$ and $C_{T}^{M}(R_{1},t)$ are independent

$\begin{array}[]{l}\begin{aligned}
E&\left[\left(C_{T}^{M}(R_{1})+C_{T}^{M}(R_{1},t)\right)^{2},R_{1}\leq
t\right]\\\ =&E\left[C_{T}^{M}(R_{1})^{2},R_{1}\leq t\right]\\\
+&E\left[C_{T}^{M}(R_{1},t)^{2},R_{1}\leq t\right]\\\ +&2\leavevmode\nobreak\
E\left[C_{T}^{M}(R_{1}),R_{1}\leq t\right]E\left[C_{T}^{M}(R_{1},t),R_{1}\leq
t\right].\end{aligned}\end{array}$

Thus,

$\begin{array}[]{l}\begin{aligned} E&\Big{[}C_{T}^{M}(R_{1})^{2},R_{1}\leq
t\Big{]}\\\ =&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(R_{1})^{2},R_{1,c}=kT\right]\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(R_{1})^{2},R_{1,p}=kT\right],\end{aligned}\end{array}$

being

$\begin{array}[]{l}\begin{aligned}
E&\left[\left(C_{T}^{M}(R_{1})\right)^{2},R_{1,c}=kT\right]\\\
=&\left(C_{c}+C_{I}(k-1)+C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\right)^{2}P_{R_{1,c}}^{M}(kT),\end{aligned}\end{array}$

and

$\begin{array}[]{l}\begin{aligned}
E\left[\left(C_{T}^{M}(R_{1})\right)^{2},R_{1,p}=kT\right]=\left(C_{c}+C_{I}(k-1)\right)^{2}P_{R_{1,p}}^{M}(kT).\end{aligned}\end{array}$

Following the same reasoning as in Appendix A,

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(R_{1}),R_{1}\leq
t\right]\\\ =&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)\right)P_{R_{1,c}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)P_{R_{1,p}}^{M}(kT).\\\
\end{aligned}\end{array}$

Since $C_{T}^{M}(R_{1},t)$ is stochastically the same as $C_{T}^{M}(t-R_{1})$,
we get

$\begin{array}[]{l}\begin{aligned} E\left[C_{T}^{M}(R_{1},t),R_{1}\leq
t\right]=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(t-kT)\right]P_{R_{1}}^{M}(kT).\end{aligned}\end{array}$

Thus,

$\begin{array}[]{l}\begin{aligned} E&\left[C_{T}^{M}(R_{1}),R_{1}\leq
t\right]E\left[C_{T}^{M}(R_{1},t),R_{1}\leq t\right]\\\ =&\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,c}}^{M}(kT)\\\
+&\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,p}}^{M}(kT).\end{aligned}\end{array}$

Then, $E\left[C_{T}^{M}(t)^{2}\right]$ verifies the following recursive
equation

$E\left[C_{T}^{M}(t)^{2}\right]=\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}E\left[C_{T}^{M}(t-kT)^{2}\right]P_{R_{1}}^{M}(kT)+H_{T}^{M}(t),$

being

$\begin{array}[]{l}\begin{aligned}
H_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)^{2}P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)+C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\right)^{2}\\\
&P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{c}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}C_{d}E\left[W_{T}^{M}((k-1)T,kT)\right]\\\
&E\left[C_{T}^{M}(t-kT)\right]P_{R_{1,c}}^{M}(kT)\\\ +&2\sum_{k=1}^{\lfloor
t/T\rfloor}\left(C_{p}+C_{I}(k-1)\right)E\left[C_{T}^{M}(t-kT)\right]\\\
&P_{R_{1,p}}^{M}(kT)\\\ +&\Big{(}\lfloor t/T\rfloor
C_{I}+C_{d}E\left[W_{T}^{M}(\lfloor t/T\rfloor T,t)\right]\Big{)}^{2}\\\
&\Bigg{(}1-\sum_{k=1}^{\lfloor t/T\rfloor}P_{R_{1}}^{M}(kT)\Bigg{)},\\\
\end{aligned}\end{array}$

and the result holds.

## Appendix C

For $t<T$, $A_{T}^{M}(t)$ is given by

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&P\left[t<\sigma_{M_{s}},\leavevmode\nobreak\
Y>t\right]+P\left[\sigma_{M_{s}}<t<\sigma_{L},\leavevmode\nobreak\
Y>t\right]\\\
=&\bar{F}_{\sigma_{M_{s}}}(t)\bar{F}_{1}(t)+\displaystyle\int_{0}^{t}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-u)I(u,t)du.\\\
\end{aligned}\end{array}$

For $t\geq T$, $A_{T}^{M}(t)$ is conditioned to the time to the first
replacement

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\displaystyle\sum_{j=0}^{\infty}\mathbf{1}_{\\{R_{j}\leq
t<R_{j+1}\\}}\\\ &\Big{[}P\left[O(t)<L,\leavevmode\nobreak\
Y>(t-R_{j}),\leavevmode\nobreak\ R_{1}\leq t\right]\\\
+&P\left[O(t)<L,\leavevmode\nobreak\ Y>(t-R_{j}),\leavevmode\nobreak\
R_{1}>t\right]\Big{]}.\end{aligned}\end{array}$

If $R_{1}>t$

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\Big{[}P\left[t<\sigma_{M},\leavevmode\nobreak\ Y>t\right]\\\
&+P\left[\lfloor t/T\rfloor
T<\sigma_{M}<\sigma_{M_{s}}<t<\sigma_{L},\leavevmode\nobreak\ Y>t\right]\\\
&+P\left[\lfloor t/T\rfloor T<\sigma_{M}<t<\sigma_{M_{s}},\leavevmode\nobreak\
Y>t\right]\Big{]}\mathbf{1}_{\left\\{M\leq M_{s}\right\\}}\\\
+&\Big{[}P\left[t<\sigma_{M_{s}},\leavevmode\nobreak\ Y>t\right]\\\
&+P\left[\sigma_{M_{s}}<\lfloor t/T\rfloor
T<\sigma_{M}<t<\sigma_{L},\leavevmode\nobreak\ Y>t\right]\\\
&+P\left[\sigma_{M_{s}}<\lfloor t/T\rfloor T<t<\sigma_{M},\leavevmode\nobreak\
Y>t\right]\\\ &+P\left[\lfloor t/T\rfloor
T<\sigma_{M_{s}}<t<\sigma_{L},\leavevmode\nobreak\
Y>t\right]\Big{]}\mathbf{1}_{\left\\{M>M_{s}\right\\}}.\end{aligned}\end{array}$

That is

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\Big{[}\bar{F}_{\sigma_{M}}(t)\bar{F}_{1}(t)\\\
&+\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}}(u)\int_{u}^{t}f_{\sigma_{M_{s}}-\sigma_{M}}(v-u)\\\
&\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-v)I(v,t)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\ &+\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}}(u)\bar{F}_{\sigma_{M_{s}}-\sigma_{M}}(t-u)\bar{F}_{1}(t)\leavevmode\nobreak\
du\Big{]}\mathbf{1}_{\left\\{M\leq M_{s}\right\\}}\\\
+&\Big{[}\bar{F}_{\sigma_{M_{s}}}(t)\bar{F}_{1}(t)\\\
&+\displaystyle\int_{0}^{\lfloor t/T\rfloor
T}f_{\sigma_{M_{s}}}(u)\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M}-\sigma_{M_{s}}}(v-u)\\\
&\bar{F}_{\sigma_{L}-\sigma_{M}}(t-v)I(u,t)\leavevmode\nobreak\
dv\leavevmode\nobreak\ du\\\ &+\displaystyle\int_{0}^{\lfloor t/T\rfloor
T}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{M}-\sigma_{M_{s}}}(t-u)I(u,t)\leavevmode\nobreak\
du\\\ &+\displaystyle\int_{\lfloor t/T\rfloor
T}^{t}f_{\sigma_{M_{s}}}(u)\bar{F}_{\sigma_{L}-\sigma_{M_{s}}}(t-u)I(u,t)\leavevmode\nobreak\
du\Big{]}\mathbf{1}_{\left\\{M>M_{s}\right\\}}\\\
&=J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}}.\end{aligned}\end{array}$

If $R_{1}\leq t$,

$\begin{array}[]{l}\begin{aligned}
\displaystyle\sum_{j=0}^{\infty}&\mathbf{1}_{\\{R_{j}\leq
t<R_{j+1}\\}}P\left[O(t)<L,\leavevmode\nobreak\
Y>(t-R_{j}),\leavevmode\nobreak\ R_{1}\leq t\right]\\\
=\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}&P_{R_{1}}^{M}(kT)\Bigg{[}\sum_{j=0}^{\infty}\mathbf{1}_{\\{R_{j}\leq
t<R_{j+1}\\}}\\\ &P\left[O(t-kT)<L,\leavevmode\nobreak\ Y>(t-kT-
R_{j})\right]\Bigg{]}\\\ =\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}&A_{T}^{M}(t-kT)P_{R_{1}}^{M}(kT).\end{aligned}\end{array}$

Then, for $t\geq T$, $A_{T}^{M}(t)$ verifies the following recursive equation

$\begin{array}[]{l}\begin{aligned}
A_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}A_{T}^{M}(t-kT)P_{R_{1}}^{M}(kT)\\\
+&J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$

and the result holds.

## Appendix D

For $t<T$, thee is no maintenance action on $[0,t]$, hence $R_{T}^{M}(t)$ is
equal to $A_{T}^{M}(t)$.

For $t\geq T$, $R_{T}^{M}(t)$ is conditioned to the time of the first
replacement

$\begin{array}[]{l}\begin{aligned}
R_{T}^{M}(t)=&P\left[O(u)<L,\leavevmode\nobreak\ \forall
u\in(0,t],\leavevmode\nobreak\ N_{s}(0,t)=0,\leavevmode\nobreak\ R_{1}\leq
t\right]\\\ +&P\left[O(u)<L,\leavevmode\nobreak\ \forall
u\in(0,t],\leavevmode\nobreak\ N_{s}(0,t)=0,\leavevmode\nobreak\
R_{1}>t\right].\end{aligned}\end{array}$

If $R_{1}>t$

$\displaystyle R_{T}^{M}(t)$ $\displaystyle=$ $\displaystyle A_{T}^{M}(t)$
$\displaystyle=$ $\displaystyle J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}}.$

If $R_{1}\leq t$,

$\begin{array}[]{l}\begin{aligned}
R_{T}^{M}(t)=&P\left[O(u)<L,\leavevmode\nobreak\ \forall
u\in(0,t],N_{s}(0,t)=0,\leavevmode\nobreak\ R_{1}\leq t\right]\\\
=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1,p}}^{M}(kT)P\left[O(u-kT)<L,\right.\\\ &\left.\forall
u\in(0,t-kT],N_{s}(0,t-kT)=0\right]\\\ =&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}R_{T}^{M}(t-kT)P_{R_{1,p}}^{M}(kT).\end{aligned}\end{array}$

Then, for $t\geq T$, $R_{T}^{M}(t)$ verifies the following recursive equation

$\begin{array}[]{l}\begin{aligned}
R_{T}^{M}(t)=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}R_{T}^{M}(t-kT)P_{R_{1,p}}^{M}(kT)\\\
+&J_{T,1}^{M}(t)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$

and the result holds.

## Appendix E

For $(t+s)<T$, there is no maintenance action on $[0,t+s]$, hence
$IR_{T}^{M}(t,t+s)$ is equal to $R_{T}^{M}(t+s)$. For $t+s\geq T$,
$IR_{T}^{M}(t,t+s)$ is conditioned to the time of the first replacement

$\begin{array}[]{l}\begin{aligned} IR_{T}^{M}(t,t+s)=&P\left[O(u)<L,\forall
u\in(t,t+s],\leavevmode\nobreak\ \right.\\\
&\left.N_{s}(t,t+s)=0,\leavevmode\nobreak\ R_{1}\leq t\right]\\\
+&P\left[O(u)<L,\leavevmode\nobreak\ \forall u\in(t,t+s],\right.\\\
&\left.N_{s}(t,t+s)=0,\leavevmode\nobreak\ t<R_{1}<t+s\right]\\\
+&P\left[O(u)<L,\leavevmode\nobreak\ \forall u\in(t,t+s],\right.\\\
&\left.N_{s}(t,t+s)=0,\leavevmode\nobreak\ R_{1}\geq
t+s\right]\Big{]}.\end{aligned}\end{array}$

If $R_{1}\geq t+s$

$\begin{array}[]{l}\begin{aligned} IR_{T}^{M}(t,t+s)=&A_{T}^{M}(t+s)\\\
=&J_{T,1}^{M}(t+s)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t+s)\mathbf{1}_{\left\\{M>M_{s}\right\\}}.\end{aligned}\end{array}$

If $t<R_{1}<t+s$

$\begin{array}[]{l}\begin{aligned}
IR_{T}^{M}(t,t+s)=&P\left[O(u)<L,\leavevmode\nobreak\ \forall
u\in(t,t+s],\right.\\\ &\left.N_{s}(t,t+s)=0,\leavevmode\nobreak\
t<R_{1}<t+s\right]\\\ =&\displaystyle\sum_{k=\lfloor
t/T\rfloor+1}^{\lfloor(t+s)/T\rfloor}P_{R_{1,p}}^{M}(kT)P\left[O(u-kT)<L,\right.\\\
&\left.\forall u\in(0,t+s-kT],N_{s}(0,t+s-kT)=0\right]\\\
=&\displaystyle\sum_{k=\lfloor
t/T\rfloor+1}^{\lfloor(t+s)/T\rfloor}R_{T}^{M}(t+s-kT)P_{R_{1,p}}^{M}(kT).\end{aligned}\end{array}$

If $R_{1}\leq t$

$\begin{array}[]{l}\begin{aligned} IR_{T}^{M}(t,t+s)=&P\left[O(u)<L,\forall
u\in(t,t+s],\right.\\\ &\left.\leavevmode\nobreak\
N_{s}(t,t+s)=0,\leavevmode\nobreak\ R_{1}\leq t\right]\\\
=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}P_{R_{1}}^{M}(kT)P\left[O(u-kT)<L,\right.\\\ &\left.\forall
u\in(t-kT,t+s-kT],\right.\\\ &\left.N_{s}(t-kT,t+s-kT)=0\right]\\\
=&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}IR_{T}^{M}(t-kT,t+s-kT)P_{R_{1}}^{M}(kT).\end{aligned}\end{array}$

Then, for $t+s\geq T$, $IR_{T}^{M}(t,t+s)$ verifies the following recursive
equation

$\begin{array}[]{l}\begin{aligned}
IR_{T}^{M}(t,t+s)=&\displaystyle\sum_{k=\lfloor
t/T\rfloor+1}^{\lfloor(t+s)/T\rfloor}R_{T}^{M}(t+s-kT)P_{R_{1,p}}^{M}(kT)\\\
+&\displaystyle\sum_{k=1}^{\lfloor
t/T\rfloor}IR_{T}^{M}(t-kT,t+s-kT)P_{R_{1}}^{M}(kT)\\\
+&J_{T,1}^{M}(t+s)\mathbf{1}_{\left\\{M\leq
M_{s}\right\\}}+J_{T,2}^{M}(t+s)\mathbf{1}_{\left\\{M>M_{s}\right\\}},\end{aligned}\end{array}$

and the result holds.

## Acknowledgements

This research was supported by Ministerio de Economía y Competitividad, Spain
(Project MTM2015-63978-P), Gobierno de Extremadura, Spain (Project GR15106),
and European Union (European Regional Development Funds). Funding for a PhD
grant comes from Fundación Valhondo (Spain).

## References

## References

  * [1] J. van Noortwijk, A survey of the application of gamma processes in maintenance, Reliability Engineering & System Safety 94(1) (2009) 2–21.
  * [2] P. A. P. Moran, A probability theory of dams and storage systems, Aust. Jour. App. Sci. 5 (2) (1954) 116–124.
  * [3] M. Abdel-Hameed, A gamma wear process, Reliability, IEEE Transactions on Reliability 24(2) (1975) 152–153.
  * [4] Z.H. Ye, L.C. Tang, H.Y. Xu A distribution-based systems reliability model under extreme shocks and natural degradation, Reliability, IEEE Transactions on Reliability 60(1) (2011) 246–256.
  * [5] P. Sung Ho, K. Jae Hoon, Lifetime estimation of LED lamp using gamma process model, Microelectronics Reliability 57 (2016) 71–78.
  * [6] A. J. Lemoine, M. L. Wenocur, On failure modeling, Naval Research Logistics 32(3) (1985) 497–508.
  * [7] N. D. Singpurwalla, Survival in Dynamic Environments, Statistical Science 10(1) (1995) 86–103.
  * [8] I. T. Castro, N. C. Caballé, C. J. Pérez, A condition-based maintenance for a system subject to multiple degradation processes and external shocks, International Journal of Systems Science 46(9) (2015) 1692–1704.
  * [9] K. T. Huynh, I. T. Castro, A. Barros, C. Bérenguer, Modeling age-based maintenance strategies with minimal repairs for systems subject to competing failure modes due to degradation and shocks, European Journal of Operational Research 218(1) (2012) 140–151.
  * [10] N. C. Caballé, I. T. Castro, C. J. Pérez, J. M. Lanza-Gutiérrez, A condition-based maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes, Reliability Engineering & System Safety 134 (2015) 98–109.
  * [11] J. H. Cha, C. Sangüesa, I. T. Castro, Maintenance policy for a system with stochastically dependent failure modes with shock-accumulation effect, IEEE Transactions on Reliability PP (99) (2016) 1–14.
  * [12] R. Ahmad, S. Kamaruddin, An overview of time-based and condition-based maintenance in industrial application, Computers $\&$ Industrial Engineering, 63 (2012) 135–149.
  * [13] W. Li, H. Pham, A condition-based inspection-maintenance model based on geometric sequences for systems with a degradation process and random shocks, Life Cycle Reliability and Safety Engineering 1(1) (2012) 26–34.
  * [14] S. M. Asadzadeh, A. Azadeh, An integrated systemic model for optimization of condition-based maintenance with human error, Reliability Engineering & System Safety 124 (2014) 117–131.
  * [15] X. Zhao, M. Fouladirad, C. Bérenguer, L. Bordes, Condition-based inspection/replacement policies for non-monotone deteriorating systems with environmental covariates, Reliability Engineering & System Safety 95 (8) (2010) 921–934.
  * [16] W. Li, H. Pham, A condition-based maintenance model for periodically inspected systems subjected to competing failure processes, International Journal of System Assurance Engineering and Management 2 (3) (2011) 226–233.
  * [17] M. Fouladirad, A. Grall, Condition-based maintenance for a system subject to a non-homogeneous wear process with a wear rate transition, Reliability Engineering & System Safety 96(6) (2011) 611–618.
  * [18] R. Ahmadi, A new approach to modeling condition-based maintenance for stochastically deteriorating systems, International Journal of Reliability, Quality and Safety Engineering 21(5) (2014) 1450024.
  * [19] J. Wu, T. Shen, M. Xie, H.Z. Huang Analysis of maintenance policies for finite life-cycle multi-state systems, Computers $\&$ Industrial Engineering 59(4) (2010) 638–646.
  * [20] S. Taghipour, D. Banjevic, A. K. Jardine, Periodic inspection optimization model for a complex repairable system, Reliability Engineering & System Safety 95(9) (2010) 944–952.
  * [21] T. Cheng, M. D. Pandey, J. van der Weide, The probability distribution of maintenance cost of a system affected by the gamma process of degradation: Finite time solution, Reliability Engineering & System Safety 108 (2012) 65–76.
  * [22] M. Pandey, T. Cheng, J. van der Weide, Finite-time maintenance cost analysis of engineering systems affected by stochastic degradation, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability 225(2) (2011) 241–250.
  * [23] N. Limnios, Reliability measures of semi-markov systems with general state space, Methodology and Computing in Applied Probability 14 (4) (2012) 895–917.
  * [24] I. T. Castro, S. Mercier, 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) 1–22doi:10.1177/1748006X16641789.
  * [25] K. T. Huynh, A. Barros, C. Bérenguer, I. T. Castro, A periodic inspection and replacement policy for systems subject to competing failure modes due to degradation and traumatic events, Reliability Engineering & System Safety 96(4) (2011) 497–508.
  * [26] J. Bertoin, Lévy Processes, Cambridge Tracts in Mathematics, Cambridge University Press, 1998.