# Imperfect and worse than old maintenances for a gamma degradation process

Franck Corset, Mitra Fouladirad, Christian Paroissin

► To cite this version:

Franck Corset, Mitra Fouladirad, Christian Paroissin. Imperfect and worse than old maintenances for a gamma degradation process. Applied Stochastic Models in Business and Industry, 2024, ENBIS 2022, 40 (3), pp.620-639. 10.1002/asmb.2849. hal-04462980

HAL Id: hal-04462980

https://cnrs.hal.science/hal-04462980v1

Submitted on 16 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Imperfect and worse than old maintenances for a gamma degradation process

Franck Corset<fnref n="1" /> | Mitra Fouladirad<fnref n="2" /> | Christian Paroissin<fnref n="3" />

**Correspondence**

Mitra Fouladirad, Aix Marseille
Université, M2P2, UMR CNRS 7340,
Cenrale Med, Marseille, France.
Email: mitra.fouladirad@univ-amu.fr

**Funding information**

French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université - A*MIDEX, Grant/Award Number: AMX-19-IET-010

**Abstract**

This article considers a condition-based maintenance for a system subject to deterioration. The deterioration is modeled by a non-homogeneous gamma process, more precisely the gamma process and the preventive maintenance are imperfect or worse than old. The corrective maintenance actions are as good as new. The maintenance efficiency or non-efficiency parameters as well as the deterioration parameters are considered to be unknown. The monitoring data under consideration give indirect information on the maintenance parameters. Therefore, an expected maximum algorithm is applied for parameter estimation.

**KEYWORDS**

ARD effect, condition-based maintenance, degradation, EM algorithm, gamma process, imperfect maintenance, statistical inference, worse than old maintenance

# 1 | INTRODUCTION

Statistical inference in presence of monitoring and maintenance data is an important issue which is drawing lots of attention lately. In such a case, maintenance actions are not only performed at failure times or a pre-fixed dates, but can be also defined accordingly to the degradation level of the system. In the case of condition-based maintenance (CBM) policy, maintenance operations are triggered if the monitored health condition indicator or deterioration level exceeds a critical threshold, refer to Reference 1. Hence, for maintenance planning and optimization, it is essential to be able to predict the deterioration behavior. It follows that degradation modeling is becoming an important issue. Stochastic models such as gamma process, Wiener process, inverse gamma process and inverse-Gaussian process are the major tools for degradation modeling, refer to References 2–5. For these models, it is assumed implicitly that the population under consideration is homogeneous (otherwise, one has to consider models with covariates and/or random effects—known as frailty models).

The classical age-replacement (as one of the main maintenance policies) for degradation models have been substantially improved by Finkelstein et al.,⁶ but still under the assumption of replacing the item by a new one (perfect replacement). Esposito et al.⁷ has proposed and studied numerically the performance of a gamma degradation model including both an age-replacement policy and a CBM policy. Their model integrates also various sources of individual heterogeneity, a situation that many times data reflects. However, here also, the component is replaced by a new one after a maintenance action.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Authors. *Applied Stochastic Models in Business and Industry* published by John Wiley & Sons Ltd.

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, LJK, Grenoble, France</footnote>

<footnote n="2">²Aix Marseille Université, M2P2, UMR CNRS 7340, Cenrale Med, Marseille, France</footnote>

<footnote n="3">³Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France</footnote>
Degradation models integrating a CBM policies have been studied in the literature.⁸ For instance, Hao et al.⁹ have considered the situation of a CBM with imperfect inspections for continuous degradation processes. In their article, they have assumed that both corrective and preventive maintenance (PM) put the system to an as-good-as new (AGAN) state. More recently, Zhao et al.¹⁰ have considered a degradation model (time-transformed Wiener process) for warranted products with imperfect maintenance, but these imperfect maintenance are deterministic in order to make the system re-starting from a required objective degradation level. Ma et al.¹¹ have considered a multi-phase Wiener process to model degradation and maintenance actions: the degradation process enters in a new phase after a maintenance action and it leads to new parameters for drift and volatility. Thus, the effect of the imperfect maintenance is both described through a reduction of the degradation level by a certain factor and through new values of the parameters for the degradation process. Beside, they have assumed that all these factors and parameters are random. As in many papers dealing with imperfect maintenance for degradation model, they consider that maintenance actions are performed at fixed and known instants. The case of preventive and opportunistic maintenance has been studied by Wang et al.¹² who have recently considered the case of a heterogeneous balanced system composed of multiple interchangeable components, the system being subject to three competing failures. They have developed a semi-Markov decision process (SMDP) to determine the optimal policy by minimizing the average maintenance cost. The case of a series system with two non-identical units have been also considered by Wang et al.¹³ under both a condition-based and an age-based replacement policy, but also with an inventory policy. Cheng and Zhao¹⁴ have recently analysed the problem of an optimal maintenance for dependent two-component degrading systems subject to imperfect repair.

In the framework of these models, parameter estimation in presence of monitoring or historical data have been addressed in literature, refer to References 5, 15–19. Degradation model calibration and parameter estimation in the case of maintenance operation has been widely considered, for instance refer to References 15, 16, and 20. As a related issue, Giorgio and Pulcini²¹ have recently studied the effect of model mis-specification of the bounded transformed gamma process on maintenance optimization. However, the joint maintenance and deterioration parameter estimation has been scarcely discussed, for instance refer to References 18, 22, and 23. In this context, the available information is very rich in comparison to the lifetime observations but the parameters estimation is also more challenging due to the complexity of the degradation phenomenon under consideration. Under the assumption of perfect maintenance actions, the maintenance parameter estimations can be derived easier than in the case of imperfect maintenance where usually the imperfect maintenance actions restore the system to an unknown level, refer to References 18, 22–24. Kamranfar et al.²³ consider truncated gamma distributed reduction parameter where the parameters are to be estimated. Corset et al.²⁴ consider an unknown reduction parameter in the framework of an ARD$_{∞}$ model and the joint estimation of the maintenance and model parameters is addressed where only the degradation level before maintenance is observed.

In this article, similarly to References 23 and 24, the system deterioration is modeled by a non-homogeneous gamma process.²,³ The statistical inference of such a process has been studied from several point of view and for several configurations.²,³,²⁴⁻²⁸ Moreover, similarly to Reference 24, both perfect and imperfect maintenance through an ARD$_{∞}$ model are considered which means after maintenance, the system is not necessarily restored to as good as new. The model considered here is quite different to the one studied by Salles et al.²² in which the authors rather focus on semi-parametric estimation for a non-homogeneous gamma process under an ARD$_{1}$ model for the maintenance effect. In addition, authors considered that a maintenance action is performed at fixed and known instants whatever the current degradation level. Regarding the different possible model for the maintenance effect, one could read the interesting work by Mercier and Castro.²⁹ In this article, authors provide stochastic comparisons for different maintenance policies and for two models of the maintenance effect (namely ARD$_{1}$ and ARA$_{1}$—arithmetic reduction of age).

The originality of our article in comparison to other ones,²⁴ resides in the consideration of “worse than old maintenance” where the deterioration level after a PM can either be reduced due to an imperfect repair or can increase due to an inappropriate maintenance operation. This latter can be due to human errors, lack of experience of the maintenance operator or harsh environmental conditions when the maintenance operations are not always implemented as it is expected. For instance, one can enumerate the case when inexperienced maintenance operators by error damage the system during a maintenance operation, the maintenance of an offshore oil platform aboard of a boat during a tempest and so forth. The “worse than old maintenance” has been considered in Reference 30 by using a mixture model where good components are replaced on failure by weak units due to lack of experience or opportunity. Cha³¹ considered a worse than old situation where maintenance modeled by a generalized Pólya process removes failures but the result is not the same as a minimal repair which the authors in Reference 32 called “as worse than old.” Cha et al.³³ has considered the “worse than minimal” repair problem which occurs, according to them, due to the adverse effects of previous repairs, but also due to other factors like environmental and internal shocks. Lee and Cha³⁴ give an example of such maintenance models. For instance,
in a multi-component system or load sharing systems, after a component failure, the working environment frequently becomes more hostile and the reliability performances become worse than before and in this case, the overall state of the system after the repair of the failed component will be worse than the state it had just before failure. Pérez-Ocón et al.³⁵ considered the case of an inexperienced maintenance operator which can deteriorate the system after repair and they proposed a semi-Markov model involving geometric processes to model the maintenance model. Babykina and Couallier³⁶ proposed to model the worse than old effect by a linear extension of Yule process initially proposed by.³⁷ In this model, contrarily to virtual age models, the conditional intensity increases essentially by the number of previous repair actions than by time. The worse than old maintenance has been mentioned or considered in different frameworks or practical cases.³⁸⁻⁴⁰ The consideration of the worse and old maintenance seems to be very new in CBM modeling.

In this article, it is assumed that it is not possible to know in advance which operation will improve or deteriorate the system. If the degradation level before and after maintenance is available the maintenance efficiency estimation is straightforward. However, in practice such situations are not always affordable or technically possible. If the degradation level is only available before or after maintenance, the indirect observation of the maintenance efficiency requires the implementation of EM algorithm for parameter estimation. Moreover, in this article, it is considered that data are collected only after maintenance actions. Therefore, the maintenance parameter should be estimated by indirect measures via an EM algorithm. These algorithms have been already discussed in maintenance modeling, refer to Reference 41. In the specific case of imperfect maintenance models, this algorithm showed good performances, refer to References 42 and 43. Considering parameter estimation with indirect observations in degradation modeling, EM algorithm has been applied with satisfying estimation results.⁴⁴,⁴⁵ More particularly, the use of the EM algorithm considering a gamma process as deterioration model has been discussed in References 17 and 46. Furthermore, the use and efficiency of the EM algorithm have been highlighted in the framework of degradation modeling and CBM in the literature, refer to References 15,16,19, and 47. The main objective of this article is to propose an estimation procedure for a non-homogeneous gamma process parameters as well as the maintenance (imperfect and as worse and old) parameters. Indeed, once the parameters are known it is possible to remediate to problems, improve maintenance planning and operations and adjust the maintenance decision rule to the field data.

The remainder of the article is as follows. In Section 2, the degradation and maintenance model as well as the maintenance cost evaluation are presented. Section 3 is devoted to the statistical inference problem where the EM algorithm is presented. In Section 4, numerical results permit to analyse the implementation of the inference method and highlight its performances. The article will end with some conclusions and perspectives.

## 2 | DEGRADATION AND MAINTENANCE MODEL

### 2.1 | Gamma process as degradation model

It is considered that the system is subject to a monotonous degradation. The degradation is modeled by a gamma process, {X_t, t ≥ 0}, which is a continuous-time stochastic process satisfying the following conditions:

* X(0) = 0 a.s.

* X(t) - X(s) ~ G(a(t) - a(s), b) for all t > s ≥ 0, the gamma distribution, with the following probability density function:

$$f_{a(t)-a(s),b}(x) = \frac{b^{a(t)-a(s)}}{\Gamma(a(t)-a(s))} x^{a(t)-a(s)-1} e^{-bx} \mathbb{I}_{x \in \mathbb{R}^{+}}, \quad (1)$$

with shape parameter a(t) = αtβ, α > 0, β > 0, and scale parameter b > 0 and Γ(.) the gamma function.

* For all t₁ < t₂ < ... < tₙ, X(tᵢ₊₁) - X(tᵢ) are independent of X(tⱼ₊₁) - X(tⱼ), ∀ 1 ≤ i, j < n, i ≠ j.

### 2.2 | Maintenance model

The stochastic process of the maintained system is denoted by {Y_t, t ≥ 0}. The system is inspected periodically at τ, 2τ, ... with the unit inspection cost C_I. The set of inspection times is defined as follows {tᵢ = iτ, i ∈ N}. One shall refer to tᵢ⁻ and
$t_i^+$ as the time just before and just after the inspection time $t_i$ respectively. The failure threshold is denoted by $M$ beyond which the system does not operate satisfactorily. This threshold is not exactly a breakdown threshold but a degradation level which is considered as not acceptable for safety or security or productivity concerns.

The failure can be only detected by inspection. The preventive threshold beyond which the preventive action is triggered is denoted by $L$, $0 < L < M$. At each inspection $t_i$, with respect to the observed degradation level $Y_{t_i}$, the following actions can be triggered:

* if $Y_{t_i} \ge M$ a corrective maintenance (CM) is performed and the system is replaced by a new one (as good as new):

- $Y_{t_i^+} = 0$ with probability one,

- a unit cost $C_C$ is incurred,

- the system is unavailable from the failure until $t_i$,

- an unavailability cost $C_u$ per time unit is induced during the unavailability period.

* if $L \le Y_{t_i} < M$, then a PM is performed. We consider here two types of PM, an imperfect one and a worse than old one:

- an imperfect PM, and more precisely, an ARD$_\infty$ model, parameterized by $\rho \in [0, 1]$, that is, the degradation level just after the inspection (at time $t_i$) is reduced by a proportional quantity of the degradation level just before the inspection: $Y_{t_i^+} = (1 - \rho)Y_{t_i^-}$. This event occurs with a probability $p$ with cost $C_P$,

- a worse maintenance modeled by an ARD$_\infty$ model, parameterized by $\rho_w > 0$, that is, the degradation level just after the inspection (at time $t_i$) is equal to: $Y_{t_i^+} = (1 + \rho_w)Y_{t_i^-}$. This event occurs with a probability $1 - p$ with cost $C_P$.

* if $Y_{t_i} < L$, no action is performed: $Y_{t_i^+} = Y_{t_i^-}$.

One shall introduce $U_i = 1_{L \le Y_{t_i} < M}$ and $B_i$ the Boolean defined as follows:

$$B_i = \begin{cases} 1 & \text{if the } i\text{th maintenance is imperfect with probability } p \\ 0 & \text{if the } i\text{th maintenance is worse with probability } 1-p. \end{cases}$$

Assume that $B_i$ are independent and thus are i.i.d. and Bernoulli distributed with the parameter $p$.

The system is assumed to evolve according to i.i.d. copies of the gamma process $(X_t)_{t \ge 0}$ between the inspection at time $(i-1)\tau$ and the inspection at time $i\tau$, denoted by $(X_t^{(i)})_{t \ge 0}$. Therefore, during the time interval $[0, \tau]$, $\forall t \in [0, \tau)$, $Y_t = X_t^{(1)}$. At the first inspection (at time $\tau$), the following scenarios are possible:

* if $Y_{\tau^-} < L$, no action is performed and then $Y_{\tau} = Y_{\tau^-} = X_{\tau^-}^{(1)}$.

* if $L \le Y_{\tau^-} < M$, a PM is performed as follows:

- an imperfect PM (ARD$_\infty$) is performed with a probability $p$ and then the degradation level is reduced of $\rho X_{\tau}^{(1)}$ and thus $Y_{\tau} = (1 - \rho)X_{\tau}^{(1)}$, with $\rho > 0$,

- A worse maintenance (ARD$_\infty$) is performed with a probability $1-p$ and then the degradation level just after the inspection is $Y_{\tau} = (1 + \rho_w)X_{\tau}^{(1)}$, with $\rho_w > 0$,

* if $Y_{\tau^-} > M$ then the system is replaced by a new one.

After the first inspection without replacement, we have $\forall t \in [\tau, 2\tau)$,

$$Y_t = Y_{\tau} + X_t^{(1)} - X_{\tau}^{(1)} = (1 - \rho)^{U_1 B_1} (1 + \rho_w)^{U_1 (1 - B_1)} X_{\tau}^{(1)} + X_t^{(2)} - X_{\tau}^{(2)}.$$

At the second inspection (at time $2\tau$), the following scenarios are possible:

* if $Y_{2\tau^-} < L$, no action is performed, and

$$Y_{2\tau} = Y_{2\tau^-} = (1 - \rho)^{U_1} X_{\tau}^{(1)} + X_{2\tau}^{(2)} - X_{\tau}^{(2)},$$
* if $L \le Y_{2\tau^{-}} < M$, we consider two cases:

- an imperfect PM ($ARD_{\infty}$) is performed, with probability $p$, and then the degradation level is reduced of $\rho(X_{2\tau}^{(1)} - X_{\tau}^{(1)})$ and thus

$$Y_{2\tau} = Y_{2\tau^{-}} - \rho(X_{2\tau^{-}}^{(2)} - X_{\tau}^{(2)}) = (1 - \rho)^{U_{1}B_{1}}(1 + \rho_{w})^{U_{1}(1-B_{1})}X_{\tau}^{(1)} + (1 - \rho)(X_{2\tau^{-}}^{(2)} - X_{\tau}^{(2)}),$$

- a worse maintenance ($ARD_{\infty}$) is performed, with probability $1-p$, and then the degradation level is added of $\rho_{w}(X_{2\tau}^{(1)} - X_{\tau}^{(1)})$ and thus

$$Y_{2\tau} = Y_{2\tau^{-}} + \rho_{w}(X_{2\tau^{-}}^{(2)} - X_{\tau}^{(2)}) = (1 - \rho)^{U_{1}B_{1}}(1 + \rho_{w})^{U_{1}(1-B_{1})}X_{\tau}^{(1)} + (1 + \rho_{w})(X_{2\tau^{-}}^{(2)} - X_{\tau}^{(2)}),$$

* if $Y_{2\tau^{-}} > M$ then the system is replaced.

After $i$ inspections without replacement, by induction, one can deduce $\forall t \in [i\tau, (i+1)\tau)$

$$Y_t = Y_{i\tau} + (X_t^{(i+1)} - X_{i\tau}^{(i+1)}),$$

and if we denote $i = [t/\tau]$, then we get

$$Y_t = \sum_{j=1}^{i} (1 - \rho)^{U_{j}B_{j}} (1 + \rho_{w})^{U_{j}(1-B_{j})} (X_{j\tau}^{(j)} - X_{(j-1)\tau}^{(j)}) + X_{t}^{(i+1)} - X_{i\tau}^{(i+1)}.$$

Finally, degradation at inspection time is

$$Y_{(i+1)\tau} = Y_{i\tau} + \sum_{j=1}^{i+1} (1 - \rho)^{U_{j}B_{j}} (1 + \rho_{w})^{U_{j}(1-B_{j})} (X_{j\tau}^{(j)} - X_{(j-1)\tau}^{(j)}). \quad (2)$$

## 2.3 | Maintenance cost evaluation

As mentioned before, the unit inspection, preventive, CM, unavailability costs are denoted by $C_I$, $C_P$, $C_C$, and $C_u$ respectively. At time $t$, the total maintenance cost is given as follows:

$$C(t) = C_I \left[ \frac{t}{\tau} \right] + C_P \sum_{l=1}^{[t/\tau]} 1_{L \le Y_{l\tau^{-}} < M} + C_C \sum_{l=1}^{[t/\tau]} 1_{Y_{l\tau^{-}} \ge M} + C_u \sum_{l=1}^{[t/\tau]} (T_f - i\tau) 1_{Y_{l\tau^{-}} \ge M, i\tau \le T_f \le (i+1)\tau},$$

where $T_f$ is the failure time

$$T_f = \inf\{t > 0, Y_t \ge M\},$$

with the probability density function $f_{T_f}$.

Since after a CM the system is as good as new, applying the renewal theory, the long run average maintenance cost can be obtained on a renewal cycle where the length is defined by $T$ (the period between two corrective replacements)

$$T = \min\{t_i > 0, Y_{t_i} \ge M, \forall i \in \mathbb{N}\}.$$

However, this formulation is very difficult since a random number of imperfect maintenance is carried out in each renewal cycle.
To overcome this difficulty, one can calculate the long run average maintenance cost between two inspections if the Markov renewal property of the maintained degradation process is satisfied. Let be $[t_{i-1}^-, t_i^-]$ the Markov renewal cycle. Consider the Markov process describing the system state just before each maintenance action $(Z_i = (Y_{i-1}^-, t_i^-))_{i \in \mathbb{N}}$. Let $\pi$ be the steady-state distribution of $(Z_i = (Y_{i-1}^-, t_i^-))_{i \in \mathbb{N}}$, if it exists. The long run average maintenance cost can be derived as follows:

$$ \lim_{t \to \infty} \frac{C(t)}{t} = \frac{\mathbb{E}_{\pi}(C(\tau))}{\mathbb{E}_{\pi}(\tau)} = \frac{C_I + C_P \mathbb{P}_{\pi}(L \le Y_{\tau-} < M) + C_C \mathbb{P}_{\pi}(Y_{\tau-} \ge M) + C_u \int_0^{\tau} tf_{T_{f,x}}(dt)\pi(dz)}{\tau}, \quad (3) $$

where $\mathbb{E}_{\pi}$ and $\mathbb{P}_{\pi}$ are the expectation and the probability under the steady-state distribution $\pi$ respectively and $f_{T_{f,x}}$ is the probability density function associated to

$$ T_{f,x} = \inf\{t > 0, X_{t+h} > M | X_h = x\}, $$

with

$$ f_{T_{f,x}}(t) = \frac{\partial \mathbb{P}(T_{f,x} \le t)}{\partial t}, $$

where

$$ \mathbb{P}(T_{f,x} > t) = \mathbb{P}(X_t < M | X_0 = x) = \frac{\Gamma(a(t), (M-x)b)}{\Gamma(a(t))}, $$

With respect to the steady-state distribution $\pi$, the corrective and PM the probability during a Markov renewal cycle are as follows:

$$ \mathbb{P}_{\pi}(L \le Y_{\tau-} < M) = \int_{L}^{M} \int_{[0,M] \times \mathbb{R}^{+}} f_{a(\tau+t)-a(t),b}(dy-x)\pi(dx, dt), $$

$$ \mathbb{P}_{\pi}(Y_{\tau-} \ge M) = \int_{0}^{L} \int_{[M, +\infty[ \times \mathbb{R}^{+}} f_{a(\tau+t)-a(t),b}(dy-x)\pi(dx, dt). $$

The distribution $\pi$ is obtained by numerical iterations and the integrals are calculated by Monte Carlo simulations.

## 2.4 | Stationary distribution of the maintained system

Let be $q$ the transition density of the Markov process $(Z_i = (Y_{i-1}^-, t_i^-))_{i \in \mathbb{N}}$, for $t = t_i^-$ and $s = t_{i-1}^-$

$$ q(B, t, x, s) = \mathbb{P}(Y_t \in B | Y_s = x) = \int_B f_{a(t)-a(s),b}(y-x)1_{L<x} dy + p \int_B f_{a(t)-a(s),b}(y - (1-\rho)x)1_{L \le x < M} dy \\ + (1-p) \int_B f_{a(t)-a(s),b}(y - (1+\rho_w)x)1_{L \le x < M} dy + \int_B f_{a(\tau),b}(y)1_{x>M} dy, \quad (4) $$

where $f_{a,b}$ is the density function of gamma distribution with shape parameter $a$ and scale parameter $b$ defined in equation (1).

For the non-homogeneous case, as explained in References 48 and 49, under some specific hypotheses the process will be K-ergodic or ergodic and therefore an asymptotic regularity behavior can be considered and a distribution can be derived. However, as it is not possible to give a generalization in the non-homogeneous case, in this article only the cost calculations steps for the homogeneous case where $\beta = 1$ are discussed. In this case, the stationary distribution does not depend on time, let us denote $\pi$ the stationary distribution which is the solution of:

$$ \tilde{\pi}(\cdot) = \int_{[0,L]} q(\cdot, x) \tilde{\pi}(dx). \quad (5) $$
Therefore,

$$ \tilde{\pi}(dy) = \int_{0}^{\min(L,y)} \tilde{\pi}(dx) f_{\alpha\tau,b}(dy-x)) + p \int_{\min\left(\frac{y}{1-\rho}, L\right)}^{\min\left(\frac{y}{1-\rho}, M\right)} \tilde{\pi}(dx) f_{\alpha\tau,b}(dy - (1-\rho)x) \\ + (1-p) \int_{\min\left(\frac{y}{1+\rho_w}, L\right)}^{\min\left(\frac{y}{1+\rho_w}, M\right)} \tilde{\pi}(dx) f_{\alpha\tau,b}(dy - (1+\rho_w)x) + f_{\alpha\tau,b}(y) \int_{M}^{+\infty} \tilde{\pi}(dx). \quad (6) $$

The fixed-point iteration algorithm can be used to solve the Equation (6) by considering that $\pi(x)$ is the solution of $\pi_{\tau}(y) = g(\pi_{\tau}(y))$, where g is a continuous function. Considering the initialization $w_1(x) = x \exp(-x)$ the algorithm leads to the following approximation for all y at iteration k,

$$ w_k(y) = \int_{0}^{\min(L,y)} w_{k-1}(x) f_{\alpha\tau,b}(y-x) dx + f_{\alpha\tau,b}(y) \int_{M}^{+\infty} w_{k-1}(x) dx \\ + p \int_{\min\left(\frac{y}{1-\rho}, L\right)}^{\min\left(\frac{y}{1-\rho}, M\right)} w_{k-1}(x) f_{\alpha\tau,b}(y - (1-\rho)x) dx + (1-p) \int_{\min\left(\frac{y}{1-\rho}, L\right)}^{\min\left(\frac{y}{1+\rho_w}, M\right)} w_{k-1}(x) f_{\alpha\tau,b}(y - (1+\rho_w)x) dx. \quad (7) $$

The stationary distribution given above cannot be solve analytically and a numerical has to be used to solve the integral equations.

## 3 | STATISTICAL INFERENCE

Consider that data is available on a maintained system according to the maintenance policy defined in the previous section. It is supposed that in the available data the type of maintenance action (PM or CM) and the inspection times are known, the types of PM (imperfect PM and worse PM), the maintenance efficiencies and the degradation parameters are unknown. This section considers the estimation of both degradation parameters $\alpha, \beta$, and $b$ and maintenance parameters $\rho$ and $\rho_w$ which correspond respectively to the imperfect and worse maintenance and $p$ the probability that a PM is imperfect. The inspection interval $\tau$ is supposed to be known. It is supposed that the system is observed just before performing a maintenance action and the degradation level after a PM is not observed. Which means, for the inference, the observations are $Y_{t_0^{-}} = 0, Y_{t_1^{-}}, Y_{t_2^{-}}, Y_{t_3^{-}}, \dots, Y_{t_n^{-}}$ with $t_i = i\tau$ for any $i \in \{0, \dots, n\}$.

### 3.1 | Likelihood function

The system is observed at times $t_i = i\tau$ for $i \in \{0, \dots, n\}$. Let be $y_0 = 0$ and $y_1, \dots, y_n$ the observed degradation levels corresponding to $Y_{t_1^{-}}, \dots, Y_{t_n^{-}}$. Recall that, after a CM, the system both degradation level and age are reset to zero. Let be $N$ the number of CM performed between 0 and $n\tau$ and denote by $t_{i_1}, \dots, t_{i_N}$ the CM times.

As it is underlined in Reference 24, the corrective replacement leads to consider a modified time scale ($s_1, \dots, s_n$) defined as follows:

* if $N = 0$ (no CM is performed on the observed time span), $s_i = t_i$ for all $i \in \{0, \dots, n\}$,

* if $N \ge 1, \forall k \in \{0, \dots, N\}, \forall i \in \{i_k, \dots, i_{k+1}\}, s_i = t_i - t_{i_k}$.

Therefore, the random variable $\Delta Y_i = Y_{t_i^{-}} - Y_{t_{i-1}^{-}}$, for any $i \in \{1, \dots, n\}$, follows a gamma distribution with shape parameter

$$ \Delta a_i = a(s_i) - a(s_{i-1}) = \alpha (s_i^\beta - s_{i-1}^\beta) = \alpha \sum_{k=0}^{N} \{(t_i - t_{i_k})^\beta - (t_{i-1} - t_{i_k})^\beta\} \mathbb{I}_{i_k \le i \le i_{k+1}}, $$
and scale parameter b. One shall denote by $f_{\Delta a_i, b}$, the probability density function of $\Delta Y_i$. Moreover, these increments are independent. It is important to highlight that only $y_i$, the realization of $Y_{t_i^-}$, is observed and $Y_{t_{i-1}^+}$ is not observed but it can be derived by $Y_{t_{i-1}^-}$ as follows:

* if $Y_{t_{i-1}^-} \le L$, no maintenance action is performed and thus, thanks to the stochastic continuity of the gamma process, we have $\mathbb{P}[Y_{t_{i-1}^+} = Y_{t_{i-1}^-}] = 1$ and hence $Y_{t_{i-1}^+} \approx Y_{t_{i-1}^-}$,

* if $L \le Y_{t_{i-1}^-} \le M$, an imperfect maintenance action is performed, which could reduce the degradation level with probability $p$ or increase the degradation level with probability $1-p$. Thus, similarly as above, we can write that $Y_{t_{i-1}^+} \approx (1-\rho)Y_{t_{i-1}^-}$ with probability $p$ and $Y_{t_{i-1}^+} \approx (1+\rho_w)Y_{t_{i-1}^-}$ with probability $1-p$,

* if $Y_{t_{i-1}^-} \ge M$, a perfect maintenance action is performed and thus the system restart as a new one, so that $Y_{t_{i-1}^+} = 0$.

Therefore, the likelihood function $L$ can be derived by the latter and the independence of increments as follows:

$$L(\alpha, \beta, b, \rho, \rho_w, p; s_0, y_0, s_1, y_1, \dots, s_n, y_n) = \prod_{i=1}^{n} f_{\Delta a_i, b}(y_i - (1 - \rho_{i-1})y_{i-1}),$$

where $\rho_{i-1}$ is the maintenance efficiency. This efficiency is considered to be a random distributed according to a mixture of two Dirac mass and a discrete random variable as follows:

$$\rho_{i-1} = \begin{cases} 0 & \text{if } y_{i-1} \le L, \\ \rho & \text{if } L \le y_{i-1} \le M \text{ with probability } p, \\ \rho_w & \text{if } L \le y_{i-1} \le M \text{ with probability } 1-p, \\ 1 & \text{if } y_{i-1} \ge M. \end{cases}$$

Denote by $\Delta_i^P = \mathbb{I}_{Y_{t_i^-} \ge M}$ the indicator of perfect maintenance. Hence,

$$Y_{t_{i-1}^-} - (1 - \Delta_i^P)(1 - \rho)^{U_i B_i} (1 + \rho_w)^{U_i (1 - B_i)} Y_{t_{i-1}^-} \sim \mathcal{G}(a(t_{i+1}) - a(t_i), b).$$

As $B_i$ are not observed, and based on the above remarks and notations, the complete likelihood function $L$ can be written as follows:

$$L(\alpha, \beta, b, \rho, \rho_w, p; s_0, y_0, s_1, y_1, \dots, s_n, y_n) = \prod_{i=1}^{n} f_{\Delta a_i, b}(y_{t_i^-} - (1 - \Delta_{i-1}^P)(1 - \rho)^{u_{i-1} b_{i-1}} (1 + \rho_w)^{u_{i-1} (1 - b_{i-1})} y_{t_{i-1}^-}),$$

where $Y_{t_0^-} = \Delta_0^P = u_0 = b_0 = 0$ and $f_{\Delta a_i, b}$ is the probability density function of the gamma distribution with shape parameter $\Delta a_i$ and scale parameter $b$.

In order to estimate the six parameters, denoted $\theta = (\alpha, \beta, b, \rho, \rho_w, p)$, we will write the complete likelihood function associated to the observations from sampling scheme considered here. As $(b_1, \dots, b_{n-1})$ are not observed, we apply the expectation-maximization algorithm, described in the next section.

## 3.2 | EM algorithm

For $i = 2, \dots, n$, we denote $\tilde{p}_i = \mathbb{P}(B_{i-1} = 1|Y_i, H_{i-1}, \theta^{(k)})$, the probability that the $(i-1)$th PM is efficient. We denote $H_{i-1} = \{Y_1, \dots, Y_{i-1}, U_1, \dots, U_{i-1}, B_1, \dots, B_{i-2}\}$, the past of the process at time $t_{i-1}$. We can write down the expectation of conditional log-likelihood (with $u_0 = 0$ and $y_0 = 0$), denoted $\ell$, as follows:

$$\ell = \sum_{i=1}^{n} \Delta a_i \log(b) - \log(\Gamma(\Delta a_i)) + \sum_{\substack{i=1 \\ u_{i-1}=0}}^{n} (\Delta a_i - 1) \log(y_i - (1 - \Delta_{i-1}^p)y_{i-1}) - b \sum_{\substack{i=1 \\ u_{i-1}=0}}^{n} (y_i - (1 - \Delta_{i-1}^p)y_{i-1})$$
$$ 
\begin{aligned}
&+ \sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} \tilde{p}_i (\Delta a_i - 1) \log(y_i - (1 - \Delta_{i-1}^P)(1 - \rho)y_{i-1}) + \sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} (1 - \tilde{p}_i) (\Delta a_i - 1) \log(y_i - (1 - \Delta_{i-1}^P)(1 + \rho_w)y_{i-1}) \\
&- b \sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} \tilde{p}_i (y_i - (1 - \Delta_{i-1}^P)(1 - \rho)y_{i-1}) - b \sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} (1 - \tilde{p}_i) (y_i - (1 - \Delta_{i-1}^P)(1 + \rho_w)y_{i-1}),
\end{aligned}
$$

where $\Gamma(.)$ is the gamma function.

We propose the following adaptation of EM algorithm:

* We maximize $\ell$ with respect to $(\alpha, \beta, b)$ by keeping only the observations $y_i$, such that $u_{i-1} = 0$ (the first line of equation). The solution is denoted $(\hat{\alpha}, \hat{\beta}, \hat{b})$.

* Thus, we maximize $\ell$ with respect to $\rho$ by keeping only the observations $y_i$ such that $u_{i-1} = 1$ and $y_i < y_{i-1}$, which ensures that the previous maintenance was effective. The solution is denoted $\hat{\rho}$.

* Finally, we apply the EM algorithm for $\rho_w$ and $p$.

At each step of the EM algorithm, we need to compute the $\tilde{p}_i$ for $i = 2, \dots, n$ such that $u_{i-1} = 1$

$$ 
\begin{aligned}
\tilde{p}_i &= P(B_{i-1} = 1 | Y_i, \mathcal{H}_{i-1}, \theta^{(k)}) \\
&= \frac{P(Y_i | B_{i-1} = 1, \mathcal{H}_{i-1}, \theta^{(k)}) P(B_{i-1} = 1 | \mathcal{H}_{i-1}, \theta^{(k)})}{P(Y_i | B_{i-1} = 1, \mathcal{H}_{i-1}, \theta^{(k)}) P(B_{i-1} = 1 | \mathcal{H}_{i-1}, \theta^{(k)}) + P(Y_i | B_{i-1} = 0, \mathcal{H}_{i-1}, \theta^{(k)}) P(B_{i-1} = 0 | \mathcal{H}_{i-1}, \theta^{(k)})} \\
&= \frac{P^{(k)} P(Y_i | B_{i-1} = 1, \mathcal{H}_{i-1}, \theta^{(k)})}{P^{(k)} P(Y_i | B_{i-1} = 1, \mathcal{H}_{i-1}, \theta^{(k)}) + (1 - P^{(k)}) P(Y_i | B_{i-1} = 0, \mathcal{H}_{i-1}, \theta^{(k)})} \\
&= \frac{P^{(k)} f_{\Delta a_i, b} (y_i - (1 - \Delta_{i-1}^P)(1 - \rho)y_{i-1})}{P^{(k)} f_{\Delta a_i, b} (y_i - (1 - \Delta_{i-1}^P)(1 - \rho)y_{i-1}) + (1 - P^{(k)}) f_{\Delta a_i, b} (y_i - (1 - \Delta_{i-1}^P)(1 + \rho_w)y_{i-1})},
\end{aligned}
$$

where $\theta^{(k)} = (\rho_w^{(k)}, p^{(k)})$ is the current value of parameter at step $k$ and $f_{\Delta a_i, b}$ is the probability density function of the gamma distribution:

$$ f_{a,b}(x) = \frac{b^a}{\Gamma(a)} x^{a-1} e^{-bx} \mathbb{1}_{x \ge 0}. $$

We update the estimation of parameters as follows:

$$ p^{(k+1)} = \frac{\sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} \tilde{p}_i}{\sum_{\substack{i=2 \\ u_{i-1}=1}}^{n} 1}, $$

and $\rho_w^{(k+1)}$ is the solution of the maximization of $\ell$.

Note that at the end of EM algorithm, we have some estimations of $\tilde{p}_i$, and we can identify a posteriori the worse CM for instance when $\tilde{p}_i$ is less than a fixed threshold.

## 4 | NUMERICAL IMPLEMENTATION

In this section, we implement simulations with different cases: homogeneous and non-homogeneous gamma processes (convex and concave cases), from practical cases (see the survey of van Noortwijk⁵⁰):

* Homogeneous case: corrosion of reinforcement ($\beta = 1$),

* Convex case: sulphate attack ($\beta = 2$),

* Concave case: diffusion-controlled ageing ($\beta = 0.5$), creep ($\beta = 1/8$) or scour-hole depth ($\beta = 0.4$).
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Trajectories of the system under maintenance operations: Homogeneous gamma process $a(t) = t$ (top left), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (top right) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom). The blue line is the preventive threshold and the red points are the worse than old maintenance action times.</figcaption></figure>

In addition, we consider cases where imperfect CM is little or very effective ($\rho \in \{0.2, 0.5, 0.8\}$) and worse maintenance is little or very harmful ($\rho_w \in \{0.1, 0.3, 0.5, 0.7\}$), while varying the proportion of harmful maintenance ($p \in \{0.5, 0.7, 0.9\}$) and the number of inspections ($n \in \{50, 100, 500\}$).

In Figure 1, examples of trajectory for $n = 50$ observations of a maintained system is illustrated. The blue line is the preventive threshold and the red line the corrective threshold. The worse maintenance cases are perceivable by red points when the PM does not reduce the degradation level. The red trajectory is the path of the initial degradation without AGAN maintenance (corrective). The three cases of homogeneous, and non-homogeneous with convex shape function and a non-homogeneous with concave shape function are considered, with $b = 1$. For the homogeneous case it is supposed that $a(t) = t$, for the convex non-homogeneous case it is supposed that $a(t) = 2t^{1.2}$ and for the concave non-homogeneous case it is supposed that $a(t) = 2t^{0.8}$. The maintenance parameters used for the illustration are as follows $\rho = 0.8$, $\rho_w = 0.1$,
<figure><img src="image_4.png" /><figcaption>Histogram of $\alpha$ (for $\alpha=2$, $\beta=0.8$, $b=2$, $\rho=0.8$, $\rho_w=0.1$, $p=0.9$, $n=500$)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Histogram of $\alpha$ (for $\alpha=2$, $\beta=0.8$, $b=2$, $\rho=0.8$, $\rho_w=0.7$, $p=0.7$, $n=500$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Histogram of $\alpha$ (for $\alpha=2$, $\beta=1.2$, $b=2$, $\rho=0.8$, $\rho_w=0.1$, $p=0.9$, $n=500$)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Histogram of $\alpha$ (for $\alpha=2$, $\beta=1.2$, $b=2$, $\rho=0.8$, $\rho_w=0.5$, $p=0.9$, $n=500$)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 2 Histogram for $\hat{\alpha}$ for different $\rho_w$ and $\rho_w = 0.1$ (left) $\rho_w \ge 0.5$ (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom).</figcaption></figure>

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

Histogram of b (for $\alpha$=1, $\beta$=1, $b$=1, $\rho$=0.5, $\rho_w$=0.1, $p$=0.9, $n$=500)

<figure><img src="image_11.png" /></figure>

Histogram of b (for $\alpha$=1, $\beta$=1, $b$=1, $\rho$=0.8, $\rho_w$=0.5, $p$=0.9, $n$=500)

<figure><img src="image_12.png" /></figure>

Histogram of b (for $\alpha$=2, $\beta$=0.8, $b$=2, $\rho$=0.8, $\rho_w$=0.1, $p$=0.9, $n$=500)

<figure><img src="image_13.png" /></figure>

Histogram of b (for $\alpha$=2, $\beta$=0.8, $b$=2, $\rho$=0.8, $\rho_w$=0.7, $p$=0.7, $n$=500)

<figure><img src="image_14.png" /></figure>

Histogram of b (for $\alpha$=2, $\beta$=1.2, $b$=2, $\rho$=0.5, $\rho_w$=0.1, $p$=0.9, $n$=500)

<figure><img src="image_15.png" /></figure>

Histogram of b (for $\alpha$=2, $\beta$=1.2, $b$=2, $\rho$=0.8, $\rho_w$=0.5, $p$=0.9, $n$=500)

<caption>FIGURE 3 Histogram for $\hat{b}$ for different $\rho_w$ ($\rho_w$ = 0.1 (left) $\rho_w$ ≥ 0.5 (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom)).</caption>
<figure><img src="image_16.png" /></figure>

Histogram of $\beta$ (for $\alpha=1, \beta=1, b=1, \rho=0.5, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_17.png" /></figure>

Histogram of $\beta$ (for $\alpha=1, \beta=1, b=1, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)

<figure><img src="image_18.png" /></figure>

Histogram of $\beta$ (for $\alpha=2, \beta=0.8, b=2, \rho=0.8, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_19.png" /></figure>

Histogram of $\beta$ (for $\alpha=2, \beta=0.8, b=2, \rho=0.8, \rho_w=0.7, p=0.7, n=500$)

<figure><img src="image_20.png" /></figure>

Histogram of $\beta$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.5, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_21.png" /></figure>

Histogram of $\beta$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)

<caption>**FIGURE 4** Histogram for $\hat{\beta}$ for different $\rho_w$ and $\rho_w = 0.1$ (left) $\rho_w \ge 0.5$ (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom).</caption>
<figure><img src="image_22.png" /><figcaption>Histogram of $\rho$ (for $\alpha=1, \beta=1, b=1, \rho=0.5, \rho_w=0.1, p=0.9, n=500$)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Histogram of $\rho$ (for $\alpha=1, \beta=1, b=1, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Histogram of $\rho$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.8, \rho_w=0.1, p=0.9, n=500$)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Histogram of $\rho$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)</figcaption></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 5 Histogram for $\hat{\rho}$ for different $\rho_w$ and $\rho_w = 0.1$ (left) $\rho_w \ge 0.5$ (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom).</figcaption></figure>
<figure><img src="image_28.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=1, \beta=1, b=1, \rho=0.5, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_29.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=1, \beta=1, b=1, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)

<figure><img src="image_30.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=2, \beta=0.8, b=2, \rho=0.8, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_31.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=2, \beta=0.8, b=2, \rho=0.8, \rho_w=0.7, p=0.7, n=500$)

<figure><img src="image_32.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.8, \rho_w=0.1, p=0.9, n=500$)

<figure><img src="image_33.png" /></figure>

Histogram of $\rho_w$ (for $\alpha=2, \beta=1.2, b=2, \rho=0.8, \rho_w=0.5, p=0.9, n=500$)

<caption>FIGURE 6 Histogram for $\hat{\rho}_w$ for different $\rho_w$ and $\rho_w = 0.1$ (left) $\rho_w \ge 0.5$ (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom).</caption>
<figure><img src="image_34.png" /><figcaption>FIGURE 7 Histogram for $\hat{p}$ for different $\rho_w$ and $\rho_w = 0.1$ (left) $\rho_w \ge 0.5$ (right) for different models homogeneous gamma process $a(t) = t$ (top), non-homogeneous $a(t) = 2 \cdot t^{0.8}$ (center) non-homogeneous $a(t) = 2 \cdot t^{1.2}$ (bottom).</figcaption></figure>
<table><caption>TABLE 1 The mean square error associated to the estimation of $\alpha$ considering different $n$, $\rho_w$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><th>ρ<sub>w</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td>0.19</td><td>0.07</td><td>0.01</td><td>0.49</td><td>0.20</td><td>0.02</td><td>0.5</td><td>0.32</td><td>0.04</td></tr><tr><th>0.5</th><td>0.31</td><td>0.12</td><td>0.01</td><td>0.23</td><td>0.10</td><td>0.01</td><td>0.017</td><td>0.12</td><td>0.31</td></tr></tbody></table>

<table><caption>TABLE 2 The mean square error associated to the estimation of $\beta$ considering different $n$, $\rho_w$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><th>ρ<sub>w</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td>0.004</td><td>0.003</td><td>0.0004</td><td>0.16</td><td>0.05</td><td>0.008</td><td>0.018</td><td>0.009</td><td>0.001</td></tr><tr><th>0.5</th><td>0.005</td><td>0.002</td><td>0.0003</td><td>0.01</td><td>0.007</td><td>0.001</td><td>0.01</td><td>0.009</td><td>0.001</td></tr></tbody></table>

<table><caption>TABLE 3 The mean square error associated to the estimation of $b$ considering different $n$, $\rho_w$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><th>ρ<sub>w</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td>0.13</td><td>0.04</td><td>0.007</td><td>0.16</td><td>0.05</td><td>0.008</td><td>0.1</td><td>0.04</td><td>0.007</td></tr><tr><th>0.5</th><td>0.12</td><td>0.03</td><td>0.006</td><td>0.13</td><td>0.04</td><td>0.007</td><td>0.09</td><td>0.04</td><td>0.005</td></tr></tbody></table>

<table><caption>TABLE 4 The mean square error associated to the estimation of $\rho$ considering different $n$, $\rho_w$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><th>ρ<sub>w</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td>0.002</td><td>0.001</td><td>0.0002</td><td>0.001</td><td>0.0002</td><td>0.00001</td><td>0.0007</td><td>0.0001</td><td>0</td></tr><tr><th>0.5</th><td>0.002</td><td>0.001</td><td>0.0002</td><td>0.001</td><td>0.0004</td><td>0.00002</td><td>0.003</td><td>0.0002</td><td>0</td></tr></tbody></table>

<table><caption>TABLE 5 The mean square error associated to the estimation of $\rho_w$ considering different $n$, $\rho_w$.</caption><thead><tr><th rowspan="2">n</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><th>ρ<sub>w</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td>7.01</td><td>1.46</td><td>0.02</td><td>21.7</td><td>24.6</td><td>19.2</td><td>16.08</td><td>16.21</td><td>16.36</td></tr><tr><th>0.5</th><td>0.06</td><td>0.04</td><td>0.003</td><td>0.02</td><td>0.003</td><td>0.0001</td><td>0.004</td><td>0.001</td><td>0.00001</td></tr></tbody></table>

$p = 0.9$. It can be easily noticed that for the convex non-homogeneous case, the corrective and PM actions are more often and therefore more worse than old maintenance is applied. In the concave non-homogeneous case the CM scarcely happen (no corrective action is occurred in the Figure 1 top right) and the number of preventive actions is less than the convex non-homogeneous case. In the homogeneous case corrective actions occur but the number of preventive actions are similar to the case of concave non-homogeneous case. Indeed, by varying the parameters and the observation horizon (data size $n$) the behavior changes but the three cases differ in similar proportions.

In Figures 2–7, the histogram of the parameters estimate are depicted for $n = 500$ observations. The left columns concern the estimations considering $\rho_w = 0.1$ and the right columns concern the estimations considering $\rho_w \ge 0.5$. It can
<table><caption>TABLE 6 The mean square error associated to the estimation of *p* considering different *n*, *ρ*<sub>w</sub>, and *p*.</caption><thead><tr><th rowspan="2">n<br>(ρ<sub>w</sub>, p)</th><th colspan="3">Convex</th><th colspan="3">Homogeneous</th><th colspan="3">Concave</th></tr><tr><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th><th>50</th><th>100</th><th>500</th></tr></thead><tbody><tr><td>(0.1, 0.5)</td><td>0.98</td><td>0.99</td><td>0.99</td><td>0.86</td><td>0.86</td><td>0.85</td><td>0.45</td><td>0.21</td><td>0.01</td></tr><tr><td>(0.1, 0.7)</td><td>0.18</td><td>0.18</td><td>0.18</td><td>0.15</td><td>0.15</td><td>0.15</td><td>0.01</td><td>0.06</td><td>0.0</td></tr><tr><td>(0.1, 0.9)</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.004</td><td>0.008</td><td>0.0005</td></tr><tr><td>(0.3, 0.9)</td><td>0.006</td><td>0.008</td><td>0.006</td><td>0.007</td><td>0.009</td><td>0.006</td><td>0.003</td><td>0.006</td><td>0.0005</td></tr><tr><td>(0.5, 0.9)</td><td>0.01</td><td>0.01</td><td>0.003</td><td>0.008</td><td>0.003</td><td>0.0007</td><td>0.007</td><td>0.002</td><td>0.0005</td></tr></tbody></table>

be noticed that for $ρ_w ≥ 0.5$ all the estimates have a satisfying behavior. Even if the estimations are very sensitive to $ρ_w$, for $n = 500$, the gamma process parameters estimate seems to give reasonable results. However, the estimations of $ρ$, $ρ_w$ and $p$ are very poor when $ρ_w = 0.1$ especially for the non-homogeneous-convex and homogeneous cases. A more extensive analysis of the results is proposed in Tables 1–6 where the relative mean square error of the estimates is given for different observation size $n$ and different worse than old maintenance parameter.

In sum, the following analysis can be derived from numerical implementations:

* The estimation with $n = 500$ and $ρ_w > 0.1$ is very efficient. In this case the data size is large and the effect of the worse than old maintenance can be distinguished more precisely to estimate maintenance parameters, right hand side of Tables 1–6.

* For gamma process parameters only the data size influences the estimation quality since the unknown parameters are estimated with data before reaching the maintenance threshold, refer to Figures 2–4 and Tables 1–3.

* For all the parameters, the estimation behavior is satisfactory even if the MSE increases when $n = 50$, refer to Tables 1–5.

* For the non-homogeneous concave case, the estimation is overall more efficient than the other cases. Indeed, in this case due to the slow degradation, the maintenance actions are less often and after exceeding the preventive threshold it is much easier to distinguish the maintenance effect, refer to Figures 1, 7, 6 and Tables 4–6.

* The homogeneous and non-homogeneous convex cases are very sensitive to $ρ_w$ and $p$. With a very small $ρ_w$ the estimation method cannot distinguish between the two maintenance actions and it considers nearly all the maintenance actions as imperfect maintenance. Which leads to a very poor estimation of $ρ_w$ and considering $p = 1$, refer to Figures 7 and 6.

* The estimation of $p$ depends on values of the couple $(ρ_w, p)$ more than just $ρ_w$ or $p$.

## 5 | CONCLUSION

In this article, the condition-based of a gradually deteriorating system is considered. The degradation is modeled by a gamma process and the PM are considered as imperfect or worse than old. The CM is as good as new. After deriving the long run average maintenance cost, the parameter estimation by maximum likelihood is considered. The deterioration maintenance parameters are estimated by an EM algorithm since the effect of the maintenance is unknown. The estimation results show very good performances for a large observation horizon and non-negligible worse than better maintenance. A sensitivity analysis to different parameters is carried out. The next step is to use the estimation results through the cost calculation and measure the impact of the estimation errors on the decision making and losses. Other optimization procedures can be used in the EM-algorithm and compared to the procedure under consideration.

This study opens avenue to new perspectives. First, prior information can be considered on maintenance parameters since some prior knowledge on the maintenance team, quality and equipment can be available. In this framework, Bayesian estimation can be carried out to estimate the maintenance parameters. In this context, dealing efficiently with hyper-parameters will be a challenging study. Afterward, the impact of the estimation uncertainty on the optimal maintenance average cost and parameters could be studied. This latter is essential for a risk-informed decision making. Another possible extension will be to consider that only the deterioration level after inspection is available (saved as historical data). In this case, a different strategy should be derived for parameter estimation.
## ACKNOWLEDGMENTS

This work received support from the French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université - A*MIDEX (AMX-19-IET-010).

## DATA AVAILABILITY STATEMENT

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

## ORCID

Christian Paroissin https://orcid.org/0000-0001-7289-3491

## REFERENCES

1. Christer AH, Wang W. A simple condition monitoring model for a direct monitoring process. *Eur J Oper Res.* 1995;82:258-269.

2. van Noortwijk JM, van Gelder PHAJM. Optimal maintenance decisions for berm breakwaters. *Struct Saf.* 1996;18(4):293-309.

3. Kahle W, Mercier S, Paroissin C. Degradation Processes in Reliability. John Wiley & Sons; 2016.

4. Doksum K, Hóyland A. Models for variable-stress accelerated life testing experiments based on wiener processes and the inverse Gaussian distribution. *Dent Tech.* 1992;34:74-82.

5. Wang X, Xu D. An inverse Gaussian process model for degradation data. *Dent Tech.* 2010;52(2):188-197.

6. Finkelstein M, Cha JH, Langston A. Improving classical optimal age-replacement policies for degrading items. *Reliab Eng Syst Saf.* 2023;236:109303.

7. Esposito N, Mele A, Castanier B, Massimiliano G. A hybrid maintenance policy for a deteriorating unit in the presence of three forms of variability. *Reliab Eng Syst Saf.* 2023;237:109320.

8. Grall A, Dieulle L, Berenguer C, Roussignol M. Continuous-time predictive maintenance scheduling for a deteriorating system. *IEEE Trans Reliab.* 2002;51(2):141-150.

9. Hao S, Yang J, Bérenguer C. Condition-based maintenance with imperfect inspections for continuous degradation processes. *App Math Model.* 2020;86:311-334.

10. Zhao X, Bin L, Jianyu X, Wang XL. Imperfect maintenance policies for warranted products under stochastic performance degradation. *Eur J Oper Res.* 2023;308:150-165.

11. Ma J, Cai L, Liao G, Yin H, Si X, Zhang P. A multi-phase wiener process-based degradation model with imperfect maintenance activities. *Reliab Eng Syst Saf.* 2023;232:109075.

12. Wang J, Liu H, Lin T. Optimal rearrangement and preventive maintenance policies for heterogeneous balanced systems with three failure modes. *Reliab Eng Syst Saf.* 2023;238:109429.

13. Wang J, Qiu Q, Wang H. Joint optimization of condition-based and age-based replacement policy and inventory policy for a two-unit series system. *Reliab Eng Syst Saf.* 2021;205:107251.

14. Cheng W, Zhao X. Maintenance optimization for dependent two-component degrading systems subject to imperfect repair. *Reliab Eng Syst Saf.* 2023;204:1109581.

15. Kim MJ, Makis V, Jiang R. Parameter estimation in a condition-based maintenance model. *Stat Probab Lett.* 2010;80(21-22):1633-1639.

16. Ghasemi A, Yacout S, Ouali MS. Parameter estimation methods for condition-based maintenance with indirect observations. *IEEE Trans Reliab.* 2010;59(2):426-439.

17. Ye ZS, Xie M, Tang LC, Chen N. Semiparametric estimation of gamma processes for deteriorating products. *Dent Tech.* 2014;56(4):504-513.

18. Kahle W. Imperfect repair in degradation processes: a Kijima-type approach. *Appl Stoch Models Bus Ind.* 2019;35(2):211-220.

19. Xu J, Liu B, Zhao X. Parameter estimation for load-sharing system subject to wiener degradation process using the expectation-maximization algorithm. *Qual Reliab Eng Int.* 2019;35(4):1010-1024.

20. Neves ML, Santiago LP, Maia CA. A condition-based maintenance policy and input parameters estimation for deteriorating systems under periodic inspection. *Comput Ind Eng.* 2011;61(3):503-511.

21. Giorgio M, Pulcini G. The effect of model misspecification of the bounded transformed gamma process on maintenance optimization. *Reliab Eng Syst Saf.* 2024;241:109569.

22. Salles G, Mercier S, Bordes L. Semiparametric estimate of the efficiency of imperfect maintenance actions for a gamma deteriorating system. *J Stat Plan Inference.* 2020;206:278-297.

23. Kamranfar H, Fouladirad M, Balakrishnan N. Inference for a gradually deteriorating system with imperfect maintenance. *Commun Stat Simul Comput.* 2021;52:5811-5829.

24. Corset F, Fouladirad M, Paroissin C. Imperfect condition-based maintenance for a gamma degradation process in presence of unknown parameters. *Proc Inst Mech Eng O J Risk Reliab.* 2022;237:546-561.

25. Bousquet N, Fouladirad M, Grall A, Paroissin C. Bayesian gamma processes for optimizing condition-based maintenance under uncertainty. *Appl Stoch Models Bus Ind.* 2015;31(3):360-379.

26. Pan Z, Feng J, Sun Q. Lifetime distribution and associated inference of systems with multiple degradation measurements based on gamma processes. *Eksploat i Niezawodn.* 2016;18:307-313.

27. Bordes L, Mercier S, Remy E, Dautrême E. Partially observed competing degradation processes: modeling and inference. *Appl Stoch Models Bus Ind.* 2016;32(5):677-696.
28. Balakrishnan N, Tsai CC, Lin CT. Gamma degradation models: inference and optimal design. In: Chen DG, Lio Y, Ng H, Tsai TR, eds. *Statistical Modeling for Degradation Data*. Springer; 2017:171-191.

29. Mercier S, Castro I. Stochastic comparisons of imperfect maintenance models for a gamma deteriorating system. *Eur J Oper Res.* 2019;273(1):237-248.

30. Badía FG, Berrade MD, Cha JH, Lee H. Optimal replacement policy under a general failure and repair model: minimal versus worse than old repair. *Reliab Eng Syst Saf.* 2018;180:362-372.

31. Cha JH. Characterization of the generalized Polya process and its applications. *Adv Appl Probab.* 2014;46(4):1148-1171.

32. Tanwar M, Rai RN, Bolia N. Imperfect repair modeling using Kijima type generalized renewal process. *Reliab Eng Syst Saf.* 2014;124:24-31.

33. Cha JH, Finkelstein M, Levitin G. Characterization of the generalized Polya process and its applications. *IMA J Manag Math.* 2023;34(1):101-122.

34. Lee H, Cha JH. New stochastic models for preventive maintenance and maintenance optimization. *Eur J Oper Res.* 2016;255(1):80-90.

35. Pérez-Ocón R, Torres-Castro I. A reliability semi-Markov model involving geometric processes. *Appl Stoch Models Bus Ind.* 2002;18(2):157-170.

36. Babykina G, Couallier V. Modelling recurrent events for repairable systems under worse than old assumption. In: Nikulin M, Limnios N, Balakrishnan N, Kahle W, Huber-Carol C, eds. *Advances in Degradation Modeling*. Springer; 2010:339-354.

37. Le Gat Y. Extending the yule process to model recurrent pipe failures in water supply networks. *Urban Water J.* 2014;11(8):617-630.

38. Berrade MD, Scarf PA, Cavalcante CA. Some insights into the effect of maintenance quality for a protection system. *IEEE Trans Reliab.* 2015;64(2):661-672.

39. Pham H, Wang H. Imperfect maintenance. *Eur J Oper Res.* 1996;94(3):425-438.

40. Van Horenbeek A, Buré J, Cattrysse D, Pintelon L, Vansteenwegen P. Joint maintenance and inventory optimization systems: a review. *Int J Prod Econ.* 2013;143(2):499-508.

41. Lefebvre Y. Using the phase method to model degradation and maintenance efficiency. *Int J Reliab Qual Saf Eng.* 2003;10(4):383-405.

42. Lim T. Estimating system reliability with fully masked data under Brown-Proschan imperfect repair model. *Reliab Eng Syst Saf.* 1998;59(3):277-289.

43. Si W, Yang Q. A generalized mixed effect Kijima model and application in optimal maintenance planning. *IEEE Trans Reliab.* 2016;65(3):1551-1561.

44. Ye ZS. On the conditional increments of degradation processes. *Stat Probab Lett.* 2013;83(11):2531-2536.

45. Zhai Q, Ye ZS. Degradation in common dynamic environments. *Dent Tech.* 2018;60(4):461-471.

46. Mercier S, Meier-Hirmer C, Roussignol M. Bivariate gamma wear processes for track geometry modelling, with application to intervention scheduling. *Struct Infrastruct Eng.* 2012;8(4):357-366.

47. Ng TS. An application of the EM algorithm to degradation modeling. *IEEE Trans Reliab.* 2008;57(1):2-13.

48. Madsen RW, Isaacson DL. Strongly ergodic behavior for non-stationary Markov processes. *Ann Probab.* 1973;1:329-335.

49. Liêu B, Dorel M. Sur le comportement asymptotique de processus de Markov non homogènes. *Studia Math.* 1967;3(28):253-274.

50. Noortwijk JM. A survey of the application of gamma processes in maintenance. *Reliab Eng Syst Saf.* 2009;94(1):2-21.

**How to cite this article:** Corset F, Fouladirad M, Paroissin C. Imperfect and worse than old maintenances for a gamma degradation process. *Appl Stochastic Models Bus Ind.* 2024;1-20. doi: 10.1002/asmb.2849