Reliability Engineering and System Safety 213 (2021) 107671

Contents lists available at ScienceDirect

Reliability Engineering and System Safety

journal homepage: www.elsevier.com/locate/ress

Modeling cavitation erosion using non-homogeneous gamma process

Q. Chatenet a , E. Remy b , M. Gagnon c , M. Fouladirad d ,∗, A.S. Tahan a
a
École de Technologie Supérieure, 1100 Rue Notre-Dame Ouest, H3C 1K3, Montréal, QC, Canada
b
EDF R&D, 6 Quai Watier, 78401 Chatou Cedex, France
c
Institut de Recherche d’Hydro Québec, 1800 Boulevard Lionel-Boulet, J3X 1S1, Varennes, QC, Canada
d
Université de Technologie de Troyes, 12 Rue Marie Curie, 10300 Troyes, France

ARTICLE INFO ABSTRACT

Keywords: Although hydroelectric generating units are highly reliable, being able to accurately model their degradation
Degradation modeling level represents a real asset in industrial and financial risk management. This paper presents and models a
Gamma process common degradation phenomenon observed on hydraulic Francis turbine runners: erosive cavitation. It gives
Bootstrap
an application of stochastic processes for degradation modeling framework in presence of real laboratory
Estimation uncertainty
experimental data. For degradation modeling, a non homogeneous gamma process is proposed. The model
Delta method
Erosive cavitation
calibration is explained and asymptotic confidence intervals for the model estimate are assessed. Because
Lifetime estimation of the limited size of available dataset, bootstrap techniques are also used to evaluate statistical estimation
uncertainties on the model parameters. These uncertainties on the degradation model are then propagated in
order to analyze how they impact the distribution of the system lifetime, characterized by the hitting time for
a given degradation threshold.

1. Introduction of water is below its vapor pressure at a given temperature. For high
velocity regions on the runner blade, the local low pressure can result
1.1. Industrial context in the formation of small vapor bubbles which then collapse suddenly.
If the created bubbles implode on solid surfaces, such as the runner
Component reliability represents a significant aspect of industrial blades, surfaces eventually end up with pitting erosion [5]. Cavitation
asset management for energy producers who operate generating units results in noise, vibration, and alteration of the surface geometry of the
in highly competitive market. In this context, being able to build robust blade. In the end it may lead to loss in the turbine efficiency [6]. Over
reliability models is crucial to the ultimate success of maintenance time, a Francis runner exposed to erosive cavitation will loose material
program [1]. For this purpose, engineers and analysts may rely on which requires costly periodic inspections and repairs. A relevant tech-
existing probabilistic models developed in literature [2,3]. For complex nique to quantify this degradation in laboratory consists of measuring
systems such as hydroelectric power stations (which includes: turbine the mass of the specimen and thus deduce the material loss during the
runner, generator, civil, . . . ), available data can be limited due to
experiment. In order to optimize maintenance operations and to reduce
the systems’ high reliability but also because of the effort needed to
periods of outage, there is a need for accurate degradation prediction.
collect degradation data. The resulting reliability model may thus suffer
While progress has been made during the last decade, computational
from a lack of accuracy and the assessment of statistical estimation
fluid dynamics simulations still face challenges to anticipate cavitation
uncertainties may become a real challenge in order to evaluate model
erosion. Indeed, recent studies on the prediction of cavitation erosion
performance and robustness. Moreover, having a good knowledge of
have shown satisfying qualitative results, but stress that future work is
these uncertainties gives the analyst a better representation of pos-
sible future system behavior and allows him, for example, to make needed to produce convincing quantitative assessments [7]. In particu-
risk-informed decisions regarding operation or maintenance policies. lar, numerical models can predict the blade regions affected by erosive
cavitation but are still not able to produce robust results about the
1.2. Degradation phenomenon cavitation intensity and consequently the resulting material loss [8].
Although the following methodology focuses on cavitation erosion, a
Cavitation erosion is the physical degradation process studied in similar approach may be applied to other cumulative and monotonic
this paper [4]. This phenomenon occurs when the local static pressure degradation mechanism.

∗ Corresponding author.
E-mail address: [email protected] (Q. Chatenet).

https://doi.org/10.1016/j.ress.2021.107671
Received 16 November 2020; Received in revised form 11 March 2021; Accepted 3 April 2021
Available online 15 April 2021
0951-8320/© 2021 Elsevier Ltd. All rights reserved.
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

on real industrial data and propose to construct not only confidence

Nomenclature intervals on model parameters but to propagate them to a reliability
R∗+ Set of positive real numbers performance criterion of interest (while most papers limit themselves
𝐃𝑜𝑏𝑠 Set of observed degradation data to the estimation of model uncertainties).
The paper is structured as follows. First, we provide a brief intro-
𝑥𝑖,𝑗 Degradation level of component 𝑗 at instant
duction of probabilistic degradation models and present the particular
𝑡𝑖,𝑗
case of the model we chose to conduct this study. Next, we present
𝛿𝑖,𝑗 Increment of degradation between two
the methodology used to construct confidence intervals on model pa-
successive instants
rameters: (i) theoretical asymptotic properties of maximum likelihood
𝑢 Gamma process rate parameter
estimators and (ii) bootstrap methods. Then, we apply the previous
𝑣 Gamma process shape parameter techniques on a real dataset and discuss the results. Finally, some of
𝜽 Set of gamma process parameters the challenges faced for degradation modeling are discussed, along with
𝑛𝑗 Path size suggestions for further work.
𝑚 Number of paths
𝑁 Total number of observations 2. Degradation models
 Likelihood function
𝓁 Log-likelihood function In health indicator prognostic, when only failure times are known,
𝛷 CDF of the standard Gaussian distribution black box models such as lifetime models are used to analyze the
behavior of the system. To take advantage of the available information
𝑃𝑌𝑞 Quantile of order 𝑞 of distribution of 𝑌
provided by deterioration data (collected during system’s monitoring
𝑇𝜌 First hitting time of degradation at level 𝜌
or inspection), stochastic processes are suitable candidates to describe
HGP Homogeneous gamma process the evolution over time of the health indicator. Indeed, degradation
NHGP Non-homogeneous gamma process data can be seen as a collection of random variables (increments of
BS Birnbaum-Saunders degradation), and can be modeled in terms of a stochastic process.
ML Maximum likelihood Another applicable tool would be the use of regression models also
MLE Maximum likelihood estimation known as ‘‘degradation path models’’ [10,11].
EB Efron bootstrap
PB Parametric bootstrap 2.1. Degradation path models
MBB Moving block bootstrap
Let us denote 𝑥𝑖,𝑗 , 𝑖 = 1 … 𝑛, with 0 = 𝑥0,𝑗 ≤ 𝑥1,𝑗 ≤ ⋯ ≤ 𝑥𝑛,𝑗 , the
DM Delta method
observed degradation level related to component 𝑗, 𝑗 = 1 … 𝑚, at instant
GSS Gamma sequential sampling
𝑡𝑖,𝑗 , where 0 = 𝑡0,𝑗 < 𝑡1,𝑗 < ⋯ < 𝑡𝑛,𝑗 . For regression models, the suc-
MSE Mean squared error cessive degradation levels are considered as independent observations
which means the dynamics of degradation and its cumulative property
over the time are overlooked, as well as the way degradation values are
collected and organized. Likewise, random effect models can be suitable
1.3. Objectives for degradation exhibiting unexplained heterogeneity within a similar
component population. These models are often used in presence of
To perform accurate degradation prediction, an analyst must rely on uncontrolled environment (which is not the case for the dataset studied
a realistic model which ideally integrates the physics of degradation, in this paper, which was obtained during laboratory experiments). In
the characteristics of the operating environment and the impact of the contrast, stochastic processes model the degradation path variability by
usage [9,10]. Accurate degradation modeling is made possible through considering the dynamics and the cumulative property of degradation
the use of degradation data collected by laboratory experiments, mon- using its increments between two successive instants. For these reasons,
itoring or during inspections of the system. These models allow the regression models are discarded in the remaining of this paper.
prediction of the lifetime of the system, characterized by the first
hitting time for a given degradation level threshold. In the case of 2.2. Stochastic processes
cavitation erosion, because of the natural variability of the physical
phenomenon and the large number of unknown influencing factors, As pointed out in the survey by Noortwijk, Markov stochastic pro-
we consider that the mechanism is random. For these reasons, we cesses are well suited for modeling the temporal variability of degra-
discarded deterministic models in favor of probabilistic ones. dation [12,13]. In the context of cumulative degradation, it is relevant
This paper aims to propose a stochastic model describing the ero- to consider processes belonging to the Markov process family as the
sive cavitation phenomenon encountered on industrial equipment. The future states depend only upon the current state and not on previous
modeling is based on real data obtained under controlled environment states [14] which is the case for two successive measurements of erosive
in laboratory. The main issue is that the dataset is very small and cavitation.
the major statistical methods applied for model calibration and vali- In the Markov process class, a distinction is made between pro-
dation are based on large size data. This problem is bypassed by the cesses having a finite or countable state space and infinite state space:
bootstrap methods which enrich the dataset by generating a large set cavitation erosion measurements lie in the infinite state case as the
of resampled data. The bootstrap behavior for this specific industrial degradation level is measured using the mass loss of the specimen
framework is studied and the advantages and limits of the bootstrap which can take theoretically an infinite number of values (uncount-
for this particular industrial case is discussed. Confidence intervals for able). Likewise, we can find discrete-time or continuous-time processes
parameters estimation are proposed. A risk-informed metric, evaluating as depicted by Karlin [15]. In the framework of our study, we con-
the influence of these confidence intervals on the system lifetime sider the physical phenomenon as a time-continuous process even if
distribution (estimating the first hitting time for a given degradation degradation observations are carried out at discrete times. Many useful
level) is proposed. This latter allows to compare the different statistical continuous-time processes can be considered to model degradation such
techniques based on asymptotic and non-asymptotic results. Although as the Brownian motion with drift (or Wiener process) [16,17], the
the methods we used in this paper are standard, this study is conducted compound Poisson process [18], the inverse Gaussian process [19]

2
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

and the gamma process [20]. A characteristic feature of the Wiener process with independent, non-negative increments having a gamma
process – in the context of degradation – is that the deterioration distribution with an identical scale parameter.
may alternately increase and decrease over time (non-monotonic), A positive random quantity 𝑋 has a gamma distribution, denoted
even with an increasing trend, called Brownian motion with drift 𝐺𝑎(𝑢, 𝑣), where 𝑢, 𝑣 > 0 are the rate parameter and the shape parameter,
(see review by Guan et al. [21] and application with measurement respectively, if its probability density function (PDF) is given by:
errors in Whitmore [22]). As the degradation studied in this paper 𝑢𝑣 𝑣−1
is cumulative and strictly monotonic, inverse Gaussian process and 𝑓𝑋 (𝑥|𝑢, 𝑣) = 𝐺𝑎(𝑥|𝑢, 𝑣) = 𝑥 exp(−𝑢𝑥)𝐼R∗ (𝑥), (1)
𝛤 (𝑣) +
gamma process are good candidates because of their independent non-
where 𝐼𝐴 (𝑥) = 1 for 𝑥 ∈ 𝐴 and 𝐼𝐴 (𝑥) = 0 for 𝑥 ∉ 𝐴, 𝛤 (𝑎) is
negative increments [23,24]. With both these two processes, paths can
the incomplete gamma function for 𝑎 > 0 [36] and R∗+ the set of
be thought as the accumulation over time of an infinite number of
positive real numbers. Note that another standard parameterization of
tiny increments. These features make these processes relevant to model
the gamma distribution in the literature is (𝛼, 𝛽) with 𝛽 = 1∕𝑢 and 𝛼 = 𝑣.
gradual damage accumulating over time. IG process was used by Ye
In that case, 𝛽 is called scale parameter and 𝛼 the shape parameter.
and Chen [25] because of its inherent flexibility and the possibility to
The gamma process {𝑋(𝑡), 𝑡 ≥ 0} with rate parameter 𝑢 > 0 and
incorporate covariates and random effects. Meanwhile, Ye and Chen
shape function 𝑣(𝑡) > 0 has the following properties:
[25] stress IG process is not well received by reliability engineers or
analysts compared to Wiener or gamma processes due to its unclear • 𝑋(0) = 0,
degradation physical meaning. More over, inverse Gaussian process • 𝑋(𝜏) − 𝑋(𝑡) ∼ 𝐺𝑎(𝑢, 𝑣(𝜏) − 𝑣(𝑡)) for all 𝜏 > 𝑡 > 0,
is often considered when Wiener or gamma process fail to model • 𝑋(𝑡) has independent non-overlapping increments.
degradation data. Other Markov processes can be found in a reliability
context such as diffusion processes. In order to better control mean Note that 𝑣 must be a non-decreasing, right-continuous, real-valued
and variance of a given degradation process, Deng et al. [26] used an function on R+ with 𝑣(0) = 0. The expectation and the variance of
Ornstein–Uhlenbeck process which lies in the diffusion process family. gamma process with parameters 𝑢 and 𝑣(𝑡) are the following:
This process shares some similarities with the Wiener process with drift, 𝑣(𝑡) 𝑣(𝑡)
such as the ability to model degradation having temporary fluctuations E(𝑋(𝑡)) = , V(𝑋(𝑡)) = . (2)
𝑢 𝑢2
while the overall trend is increasing. Because the studied degradation
The shape function 𝑣(𝑡) expresses how the mean degradation in-
is monotonic and cumulative, these processes were discarded to prefer
creases over time. A standard choice is 𝑣(𝑡) = 𝑐𝑡𝑏 with 𝑐, 𝑏 > 0 [3].
a monotonic process for the remaining of the paper. Since the aim
Empirical studies and engineering knowledge show the expected degra-
is to derive a calibration, estimation and prediction methodology in
dation over time is often proportional to a power function [12]. When
presence of small size data and highlight the interest of time dependent
𝑏 = 1, the gamma process is said to be homogeneous (HGP, or
average and volatility, gamma process as a monotonic process is used
stationary); for 𝑏 ≠ 1, the gamma process is non-homogeneous (NHGP,
to show the implementation of the whole methodology. Any other Lévy
or non-stationary). In their study, Park and Padgett [37] identified
subordinator or non stationary non decreasing stochastic process could
an appropriate value of 𝑏 based on the Akaike Information Criterion,
be chosen. which is a standard statistical indicator for model selection. In the
remaining of the paper, we chose to use a power function for 𝑣. Another
2.3. Model description standard shape function could have been the exponential one: 𝑣(𝑡) =
𝑒−𝑐𝑡 − 1, with 𝑐 > 0, but because of the degradation behavior studied in
In its homogeneous form, the gamma process has a linear average this paper (see Section 4), we decided to discard this shape function.
trend but the non-homogeneous gamma process can have convex or
concave average trends. This is more flexible to fit a given aver- 2.4. Maximum likelihood estimation of model parameters
age trend but it does not satisfy the stationarity property. Castanier
et al. [27] used it in the framework of condition-based maintenance. In order to model degradation phenomenon using gamma process,
In de Jonge [28], Deloux et al. [29] and Zhu et al. [30], the gamma appropriate statistical methods are required for the model parameters’
process is enriched with additional random shocks or/and covariates to estimation. Two of the most common techniques are the Maximum
describe a more complex behavior undergoing environmental changes. Likelihood method (ML) and the Method of Moments (MM). MM is
In presence of measurements, the fitting and calibration is one of the preferred when power 𝑏 ≠ 1 in 𝑣(𝑡) is known and for i.d. time
major issues and the parameters estimation is the main challenge. increments: in that case, the NHGP can be transformed into a HGP
The latter is addressed in Wang [31] and Kahle et al. [3] under which leads to simpler parameters’ estimation using MM [12]. MM
different hypotheses on data or model. Ling et al. [24] used and being explicit, it can be used as a starting point for ML method which
estimated parameters of a non-homogeneous gamma process to study may require an iterative optimization procedure. Meanwhile, when 𝑏
the evolution of the remaining useful life of light emitting diodes under needs to be determined, ML is more convenient than MM. In our case
stress conditions. In addition, they evaluated uncertainties of model study, preliminary statistical regression showed a good fit for power
parameter estimations in the case of shape function following an in- function (see degradation path in Fig. 3). As we had no information on
verse power law. When degradation data contain noise or measurement parameter 𝑏, we decided to use the ML method. Let us introduce this
errors, leading to non-strictly monotonic trend, Le Son et al. [32] estimation method.
used a Gibbs sampling technique to filter signal and then used a non- Let us suppose we observe 𝑚 independent paths following the same
homogeneous gamma process to model hidden degradation data for the stochastic process. Here, 𝑚 expresses the number of laboratory exper-
same purpose; Lu et al. [33] preferred the use of the Genz transform iments performed (it could also express 𝑚 different and independent
and quasi-Monte Carlo method. systems on which we observe one degradation path). Let 𝑛𝑗 denote the
Abdel-Hameed [20] was the first to use gamma process (GP) to number of successive observations taken during experiment 𝑗, 1 ≤ 𝑗 ≤
model degradation occurring randomly in time. Van Noortwijk [12] 𝑚. Then, the 𝑗th degradation path is observed 𝑛𝑗 + 1 times at instants
and Alaswad and Xiang [13], in their respective survey, noted the 𝑡0,𝑗 < 𝑡1,𝑗 < ⋯ < 𝑡𝑛𝑗 ,𝑗 . We adopt the following notations in the remaining
ability of gamma process to model gradual damage monotonically of the paper, for all 1 ≤ 𝑗 ≤ 𝑚 (see Fig. 1):
accumulating over time. In particular, Cholette et al. [34] and Nystad
et al. [35] used GP to model systems subject to erosion. The definition • 𝑡0,𝑗 = 𝑡0 = 0;
given by Van Noortwijk [12] describes gamma process as a stochastic • 𝑥0,𝑗 = 0, which means at 𝑡0 there is no degradation;

3
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

By differentiating 𝓁 with respect to the three parameters, the Max-

imum Likelihood Estimator (MLE) 𝜽̂ 𝑀𝐿 = (𝑢̂ 𝑀𝐿 , 𝑐̂𝑀𝐿 , 𝑏̂ 𝑀𝐿 ) of true un-
known parameters 𝜽 = (𝑢, 𝑐, 𝑏) consists in the solution of the following
set of three non-linear equations:
∑𝑚 ∑𝑛𝑗 𝑏̂ 𝑏̂ 𝑀𝐿 ∑𝑚 𝑏̂ 𝑀𝐿
𝑗=1 𝑡 𝑀𝐿 − 𝑡𝑖−1,𝑗
𝑖=1 𝑖,𝑗 𝑗=1 𝑡𝑛𝑗 ,𝑗
𝑢̂ 𝑀𝐿 = 𝑐̂𝑀𝐿 ∑ 𝑚 ∑ 𝑛𝑗 = 𝑐̂𝑀𝐿 ∑𝑚 , (5)
𝑗=1 𝛿
𝑖=1 𝑖,𝑗 𝑗=1 𝑥𝑛𝑗 ,𝑗
𝑛𝑗 [ ]{ ( [̂ ]) }
∑
𝑚 ∑
𝑏̂ 𝑏̂ 𝑀𝐿 𝑏 𝑏̂ 𝑀𝐿
𝑡𝑖,𝑗𝑀𝐿 − 𝑡𝑖−1,𝑗 𝜓 𝑐̂𝑀𝐿 𝑡𝑖,𝑗𝑀𝐿 − 𝑡𝑖−1,𝑗 − log(𝛿𝑖,𝑗 )
𝑗=1 𝑖=1
∑
𝑚
𝑏̂
= (𝑡𝑛𝑀𝐿
,𝑗 ) log(𝑢̂ 𝑀𝐿 ), (6)
𝑗
𝑗=1

∑
𝑚 ∑𝑗 𝑛 [̂ ]
𝑏 𝑏̂ 𝑀𝐿
𝑐̂𝑀𝐿 𝑡𝑖,𝑗𝑀𝐿 log(𝑡𝑖,𝑗 ) − 𝑡𝑖−1,𝑗 log(𝑡𝑖−1,𝑗 )
𝑗=1 𝑖=1
[ ( [̂ ]) ]
𝑏 𝑏̂ 𝑀𝐿
× 𝜓 𝑐̂𝑀𝐿 𝑡𝑖,𝑗𝑀𝐿 − 𝑡𝑖−1,𝑗 − log(𝛿𝑖,𝑗 ) (7)
Fig. 1. Data notation example.
𝑚 (
∑ )
𝑏̂ ( )
= 𝑐̂𝑀𝐿 𝑡𝑛𝑀𝐿,𝑗 log(𝑡𝑛𝑗 ,𝑗 ) log 𝑢̂ 𝑀𝐿 ,
𝑗
𝑗=1
• 𝛿𝑖,𝑗 = 𝑥𝑖+1,𝑗 − 𝑥𝑖,𝑗 : the increment of degradation between two
successive measurement instants 𝑡𝑖,𝑗 and 𝑡𝑖+1,𝑗 ; where function 𝜓(𝑎), the digamma function, is the derivative of the
𝑛 { } 𝑛 −1 { } logarithm of the gamma function:
• 𝐃𝑜𝑏𝑠 = ∪𝑚 ∪ 𝑗 (𝑡𝑖,𝑗 , 𝑥𝑖,𝑗 ) = ∪𝑚
𝑗=1 𝑖=0
∪𝑗
𝑗=1 𝑖=0
(𝑡𝑖,𝑗 , 𝑡𝑖+1,𝑗 , 𝛿𝑖,𝑗 ) : the set of
all observed degradation data, organized in 𝑚 independent paths 𝛤 ′ (𝑎) 𝑑
𝜓(𝑎) = = log 𝛤 (𝑎). (8)
with 𝑛𝑗 observations, 1 ≤ 𝑗 ≤ 𝑚. 𝛤 (𝑎) 𝑑𝑎

Hence, 𝑐̂𝑀𝐿 and 𝑏̂ 𝑀𝐿 can be either numerically computed by solving

In our case, since the degradation measurements come from ex-
Eqs. (6) and (7) (replacing 𝑢̂ 𝑀𝐿 by its expression given in (5) which is
periments performed in laboratory (see Section 4), under controlled
function of 𝑐̂𝑀𝐿 and 𝑏̂ 𝑀𝐿 ) or by first computing the MLE for 𝑐 and 𝑏
conditions, we consider the measurement error is negligible compared
and then computing 𝑢̂ 𝑀𝐿 using Eq. (5). Both methods are numerical and
to the natural variability of the degradation phenomenon (see [22]
imply an optimization using proper computing tools to solve Eqs. (5)–
for stochastic process considering measurement error). Note that for
(7). We chose to use the root finding method hybr which relies on
a given experiment 𝑗, 1 ≤ 𝑗 ≤ 𝑚, increments 𝛿𝑖,𝑗 , 0 ≤ 𝑖 ≤ 𝑛𝑗 − 1,
MINPACK’s hybrd and hybrj routines (modified Powell method, see
are independent but non-identically distributed (i.non-i.d) since they
scipy.optimize.root for more details).
depend on the (possibly non-uniform) width of time intervals [𝑡𝑖,𝑗 , 𝑡𝑖+1,𝑗 ]
between successive measurements. Because of this property, the classi-
cal graphical tools and goodness of fit tests are not applicable. Indeed, 2.5. Distribution functions of the component time to failure
in the non-homogeneous case, each degradation increment follows a
different gamma distribution described, in particular, by instants 𝑡𝑖,𝑗 A system is often considered to be failed as soon as its degradation
and 𝑡𝑖+1,𝑗 (not 𝛿𝑡𝑖 anymore). Other classes of classification tools could level is beyond a given threshold 𝜌. Let 𝑇𝜌 be the time at which failure
be studied but due to their originality, for their application their whole occurs which is the first hitting time of degradation level 𝜌: 𝑇𝜌 =
background should be recalled and defined which goes further than the inf {𝑡 ≥ 0, 𝑋(𝑡) ≥ 𝜌}. Since the degradation increments follow gamma
main objective of this paper. distributions (see Eq. (1)), as it is mentioned in [12], the cumulative
The ML method consists of identifying the model parameter 𝜽̂ = distribution function (CDF) of 𝑇𝜌 can be written as follows:
(𝑢,
̂ 𝑐, ̂ values which maximize the likelihood function  given observed
̂ 𝑏)
data 𝐃𝑜𝑏𝑠 . Intuitively, it selects the parameters’ values that make the F𝑇𝜌 (𝑡|𝑢, 𝑣(𝑡)) = P(𝑇𝜌 ≤ 𝑡) = P(𝑋(𝑡) ≥ 𝜌)
∞ (9)
data most probable. The gamma process likelihood function  is given 𝛤 (𝑣(𝑡), 𝜌)
= 𝑓𝑋(𝑡) (𝑥|𝑢, 𝑣(𝑡)) 𝑑𝑥 = ,
in Eq. (3) (see [12]). For convenience, we chose to maximize the log- ∫𝑥=𝜌 𝛤 (𝑣(𝑡))
likelihood function 𝓁 in Eq. (4).; because of strictly increasing behavior
where 𝛤 (⋅, ⋅) is the upper incomplete gamma function. Moreover since
of the log function, log-likelihood achieves its maximum value at the
the shape function 𝑣 considered in our study is differentiable, the PDF
same points as the function itself.
of 𝑇𝜌 becomes:
𝑛
∏
𝑚 ∏𝑗
( ) 𝑣′ (𝑡) ∞
(𝜽|𝐃𝑜𝑏𝑠 ) = (𝑢, 𝑐, 𝑏|𝐃𝑜𝑏𝑠 ) = 𝐺𝑎 𝛿𝑖,𝑗 |𝑢, 𝑣(𝑡𝑖,𝑗 ) − 𝑣(𝑡𝑖−1,𝑗 ) ∀𝑡 ≥ 0, 𝑓𝑇𝜌 (𝑡|𝑢, 𝑣(𝑡)) = {log(𝑧) − 𝜓(𝑣(𝑡))}𝑧𝑣(𝑡)−1 𝑒−𝑧 𝑑𝑧. (10)
𝑗=1 𝑖=1
𝛤 (𝑣(𝑡)) ∫𝑢𝜌
𝑚 𝑛𝑗
∏∏ 𝑐[𝑡𝑏 −𝑡𝑏 ] 𝑏 𝑏 Note that other expressions for the CDF of hitting time 𝑇𝜌 can be
𝑢 𝑖,𝑗 𝑖−1,𝑗 𝑐[𝑡𝑖,𝑗 −𝑡𝑖−1,𝑗 ]−1
= 𝑏 𝑏
𝛿𝑖,𝑗 𝑒−𝑢𝛿𝑖,𝑗 . found in [39]. Meanwhile, because of the gamma function behavior,
𝑗=1 𝑖=1 𝛤 (𝑐[𝑡𝑖,𝑗 − 𝑡𝑖−1,𝑗 ])
computing equation (10) can be tricky in practice [37] even if the
(3) expression gives the exact distribution of the first passage time to
( ) threshold 𝜌. To solve this issue, Park and Padgett [37] introduced
𝓁(𝜽|𝐃𝑜𝑏𝑠 ) = log (𝑢, 𝑐, 𝑏|𝐃𝑜𝑏𝑠 ) the Birnbaum-Saunders distribution [40] to approximate the above
∑
𝑚 ∑𝑗 𝑛 [ ( )] distribution in the case of HGP. This approximation of the hitting time
= 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] log(𝑢) − log 𝛤 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] (4) CDF, denoted F𝐵𝑆 , is simpler than the exact expression of F𝑇𝜌 . Paroissin
𝑗=1 𝑖=1
( ) and Salami [39] proved the approximation is also valid for non linear
+ 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] − 1 log(𝛿𝑖,𝑗 ) − 𝑢𝛿𝑖,𝑗 . gamma shape function 𝑣. Eq. (9) can be then approximated by the

4
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

Fig. 2. Approach to obtain uncertainties on both estimations of model parameter and associated hitting time distribution.

Fig. 3. (Left) Characteristic stages of the cumulative erosion-time curve, adapted from [38] - (Right) Cavitation erosion paths (𝑚 = 3) from laboratory.

following formula: 3.1. Asymptotic confidence intervals approximation using the MLE conver-
gence property
F𝑇𝜌 (𝑡|𝑢, 𝑣(𝑡)) = P(𝑇𝜌 ≤ 𝑡) ≃ F𝐵𝑆 (𝑡|𝑢, 𝑣(𝑡)),
( ) ∑𝑚
𝑢𝜌 − 𝑣(𝑡) As 𝑁 = 𝑗=1 𝑛𝑗 , the total number of observations, gets larger,
with F𝐵𝑆 (𝑡|𝑢, 𝑣(𝑡)) = −𝛷 √
𝑣(𝑡) (11) the log-likelihood function 𝓁 shows two interesting properties. First,
√ √ 𝓁(𝜽|𝐃𝑜𝑏𝑠 ) becomes more sharply peaked around 𝜽̂ 𝑀𝐿 . Second, 𝓁(𝜽|𝐃𝑜𝑏𝑠 )
⎡√ ⎛ ⎞⎤
𝑣(𝑡) 𝑢𝜌 ⎟⎥ becomes more symmetrical around 𝜽̂ 𝑀𝐿 . This asymptotic behavior
= 𝛷 ⎢ 𝑢𝜌 ⎜ − ,
⎢ ⎜ 𝑢𝜌 𝑣(𝑡) ⎟⎥ shows basically that as the sample size of observations 𝐃𝑜𝑏𝑠 grows, we
⎣ ⎝ ⎠⎦
are becoming more confident that 𝜽̂ 𝑀𝐿 lies close to the true parameter
where 𝛷 is the CDF of the standard Gaussian distribution.
𝜽. Thus, the curvature of the likelihood (i.e. the second derivative of
𝓁(𝜽|𝐃𝑜𝑏𝑠 ) with respect to 𝜽) is an important characteristic of 𝜽. This
3. Model parameters’ estimations uncertainties property is the basis for constructing confidence intervals for the un-
known parameter 𝜽. These assertions are supported by the asymptotic
Determining confidence intervals for model parameters is crucial to property of first and second derivative of 𝓁(𝜽|𝐃𝑜𝑏𝑠 ) given by the central
evaluate the confidence put in the lifetime prediction by the analyst. limit theorem (CLT) [41]. Hoadley [42] gives conditions under which
Even if stochastic processes, such as gamma process, can provide a good asymptotic convergence of MLE can be demonstrated in the case of
mathematical framework for degradation modeling, MLE (among other independent but not identically distributed (i.non-i.d.) observations
statistical estimations) is only based on limited information. In our case, which is the case for our study case (see Sections 2.4 and 4 ). To lighten
information comes from the cumulative degradation measurements, the notation, 𝓁(𝜽|𝐃𝑜𝑏𝑠 ) will be simply noted 𝓁(𝜽). In the context of MLE,
observed over a specific period of time with a few data points. Thus, the asymptotic normality derived by Rao [43] gives:
different degradation paths may result in different estimations of model
𝑑 ( )
parameters. To reflect statistical uncertainty and phenomenon variabil- (𝜽̂ 𝑀𝐿 − 𝜽) ⟶  0, 𝐼 −1 (𝜽) , (12)
ity, confidence interval (CI) including the true (unknown) parameter 𝜽
𝑑
should be then computed. Fig. 2 shows the approach used to obtain with 𝐼(𝜽) the expected Fisher information matrix. Also, ⟶ denotes
uncertainties on both estimations of model parameter and hitting time convergence in law when 𝑁 goes to infinity. Because computing the ex-
distribution from methods introduced in the following sections. pected Fisher information is not possible when its analytic expression is

5
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

unavailable, we used the observed Fisher information 𝐽 (𝜽), introduced consider 𝑔(𝜽) = 𝑔𝑞 (𝜽) = F−1𝐵𝑆
(𝑞|𝜽) in the DM. The resulting function
by Efron and Hinkley [44], which is given by the following expression: 𝑔 is given by the following equation:
( )−1 ( )
𝐽 (𝜽) = −∇2 𝓁(𝜽). (13) 𝑔𝑞 (𝜽) = 𝛷[ℎ(𝑃𝑇𝑞 |𝜽)] = ℎ−1 𝛷−1 (𝑞)|𝜽
𝜌

The gradient vector of the first partial derivatives of log-likelihood √ √ 1∕𝑏

⎛ ⎞
2𝑐𝜌𝑢 + 2𝑐 𝑒𝑟𝑓𝑐−1 (2𝑞)2 − 2𝑐 𝑒𝑟𝑓𝑐−1 (2𝑞) 4𝜌𝑢 + 2 𝑒𝑟𝑓𝑐−1 (2𝑞)2 ⎟
function 𝓁(𝜽) of the gamma process can be written: −1∕𝑏 ⎜
=2 ⎜ ⎟ ,
⎜ 𝑐2 ⎟
⎛ 𝜕𝓁(𝜽) ⎞ ⎝ ⎠
⎜ 𝜕𝑢 ⎟
∇𝓁(𝜽) = ⎜ 𝜕𝓁(𝜽) ⎟, (14) (19)
⎜ 𝜕𝑐 ⎟
⎜ 𝜕𝓁(𝜽) ⎟ √ (√ 𝑣(𝑡) √ 𝑢𝜌 )
⎝ 𝜕𝑏 ⎠ where ℎ(𝑡|𝜽) = 𝑢𝜌 𝑢𝜌
− 𝑣(𝑡) refers to Birnbaum-Saunders ap-
and the Hessian matrix with the second partial derivatives of log- proximation F𝐵𝑆 (𝑡|𝜽) = 𝛷 (ℎ(𝑡|𝜽)) in Eq. (11) and 𝑒𝑟𝑓𝑐−1 denotes the
likelihood function 𝓁(𝜽) is the following: inverse function of 𝑒𝑟𝑓𝑐 which is the complementary error function
given in the following expression:
2 𝜕 2 𝓁(𝜽) 𝜕 2 𝓁(𝜽) ⎞
⎛ 𝜕 𝓁(𝜽)
2 ∞
⎜ 𝜕𝑢 𝜕𝑢𝜕𝑐 𝜕𝑢𝜕𝑏 ⎟ 2
𝑒𝑟𝑓𝑐 (𝑧) = √ 𝑒𝑥𝑝(−𝑡2 )𝑑𝑡. (20)
∇2 𝓁(𝜽) = ⎜ 𝜕 𝓁(𝜽) 𝜕 2 𝓁(𝜽) ⎟
2
𝜋 ∫𝑧
𝜕𝓁(𝜽)
. (15)
⎜ 𝜕𝑐𝜕𝑢 𝜕𝑐 2 𝜕𝑐𝜕𝑏 ⎟
⎜ 𝜕2 𝓁(𝜽) 𝜕𝓁(𝜽) 2
𝜕 𝓁(𝜽) ⎟
⎝ 𝜕𝑏𝜕𝑢 𝜕𝑏𝜕𝑐 𝜕𝑏2 ⎠ 3.3. Bootstrap techniques
It is possible to compute confidence intervals on 𝜽 using the asymp-
totic normality of MLE. The asymptotic approximate symmetric two- In the previous section, confidence intervals were built upon the
sided confidence intervals at level 1 − 𝛼, 0 < 𝛼 < 1, for 𝜽 considering assumption of asymptotic convergence of MLE, which implies that
each of its components 𝜽𝑖 are: the total number of observations is large enough to be considered
∑𝑚
[ √( √( ] as infinite: 𝑁 = 𝑗=1 𝑛𝑗 → +∞. Because of the small size of data
) )
𝜽̂ 𝑀𝐿,𝑖 − 𝑠 ̂ −1 ; 𝜽̂ 𝑀𝐿,𝑖 + 𝑠
𝐽 (𝜽) 𝐽 ( ̂ −1
𝜽) , (16) studied in Section 4, the CLT conditions may not be satisfied. In
𝑖𝑖 𝑖𝑖
this section, we consider using several bootstrap techniques to enrich
( )
𝛼 dataset and analyze the asymptotic confidence bounds. Efron [47] and
where 𝑠 = 𝛷 1 − (e.g.: 1.96 for 1 − 𝛼 = 95% confidence level).
2 Efron and Tibshirani [48] were the firsts to introduce the principle of
bootstrapping. It consists in drawing in a reference sample to built a
3.2. Delta method large number of bootstrap samples with the same size as the original
data (see [49]) and in estimating the model parameters on each of
The delta method (DM), also called the method of propagation these samples. Especially, [50] showed this type of techniques is useful
of error, allows analyst to derive the variance of a transformation of to approximate parameters’ estimators distributions of complex models
one or more random variables [45]. DM states that the asymptotic from a small number of data.
behavior of functions over a random variable can be approximated
if the random variable is itself asymptotically normal. In particular, 3.3.1. Parametric bootstrap
DM can be used to derive percentiles of hitting time distribution from For parametric bootstrap (PB), we make the assumption that data
estimated parameters and their respective CI. To our knowledge, we are come from a given parametric model and the fitted gamma process and
the first to apply delta method in the framework of gamma processes. its associated parameter 𝜽̂ 𝑀𝐿 are used to generate new datasets [51,
∑𝑚
First, as described in [46], let us consider the sequence (𝑊𝑛 )𝑛∈N of 52]. Given the observations 𝐃𝑜𝑏𝑠 , with size 𝑁 = 𝑗=1 𝑛𝑗 , we first
R𝑝 -valued random vectors. Assume there exist 𝜇𝑊 ∈ R𝑝 and 𝛴 a 𝑝 × 𝑝 ̂ ̂
estimate 𝜽𝑀𝐿 (see Section 2.4). Then, 𝜽𝑀𝐿 is used to resample 𝐵 times
definite positive matrix such that: 𝑚 degradation paths of length 𝑛𝑗 , 1 ≤ 𝑗 ≤ 𝑚, of the gamma process using
√ the GSS method [53]. Hence, the resulting 𝐵 bootstrap samples are used
𝑑 ∗
←  (0, 𝛴).
𝑛(𝑊𝑛 − 𝜇𝑊 ) ←←←←←←←←←←→ (17) to obtain 𝐵 estimates of 𝜽̂ 𝑀𝐿 ∶ 𝜽̂ 𝑃 𝐵,𝑣 , 1 ≤ 𝑣 ≤ 𝐵. From these 𝐵 values of
𝑛→∞ ∗
𝜽̂ 𝑃 𝐵 , we can derive information on the estimation uncertainty on 𝜽̂ 𝑀𝐿
Consider 𝑟 real functions 𝑔1 , … , 𝑔𝑟 with continuous first partial deriva-
and thus on how 𝜽̂ 𝑀𝐿 behaves in comparison to 𝜽. We can illustrate
tives in 𝜇𝑊 , where at least one of these derivatives is non-zero. For the procedure with the following algorithm:
𝑖 ∈ {1, … , 𝑟} and for any 𝑛 ∈ N∗ , set 𝑍𝑖,𝑛 = 𝑔𝑖 (𝑊𝑛 ), 𝑍𝑛 = (𝑍1,𝑛 , … , 𝑍𝑟,𝑛 )𝑇
( )𝑇
and 𝜇𝑍 = 𝑔1 (𝜇𝑊 ), … , 𝑔𝑟 (𝜇𝑊 ) . Then the sequence (𝑍𝑛 )𝑛∈N∗ is also 1. Given the degradation observations 𝐃𝑜𝑏𝑠 , fit a non-homogeneous
asymptotically normal: gamma process. This leads to MLE 𝜽̂ 𝑀𝐿 .
√ 𝑑 2. Generate 𝑚 degradation paths of length 𝑛𝑗 , 1 ≤ 𝑗 ≤ 𝑚, following
←  (0, 𝐾𝛴𝐾 𝑇 ),
𝑛(𝑍𝑛 − 𝜇𝑍 ) ←←←←←←←←←←→ (18) the gamma process having 𝜽̂ 𝑀𝐿 as parameter on discrete time
𝑛→∞
interval [0, 𝑡𝑛𝑗 ].
𝑑𝑔𝑖 ∗
where 𝐾 is the 𝑟 × 𝑝 matrix with elements 𝑘𝑖,𝑗 = for 𝑖 ∈ {1, … , 𝑟} 3. Estimate 𝜽̂ 𝑃 𝐵 from step 2. new dataset.
𝑑𝑤𝑗
∗
and 𝑗 ∈ {1, … , 𝑝}. 4. Repeat steps 2. and 3. 𝐵 times to get 𝐵 MLE: 𝜽̂ , 1 ≤ 𝑣 ≤ 𝐵,𝑃 𝐵,𝑣
In this study, we are interested in how estimation uncertainties on 𝜽 from 𝐵 simulated datasets of 𝑁 measurements.
∗
may impact the distribution of 𝑇𝜌 , in particular its percentiles of order 5. From the 𝐵 obtained 𝜽̂ 𝑃 𝐵,𝑣 , 1 ≤ 𝑣 ≤ 𝐵, compute the empirical
𝑞, 0 < 𝑞 < 1, given by: 𝑃𝑇𝑞 = F−1 (𝑞|𝜽) ≈ F−1 (𝑞|𝜽). If we consider the ∗ 1 ∑𝐵 ̂ ∗
𝜌 𝑇𝜌 𝐵𝑆 mean 𝜽̂ 𝑃 𝐵 = 𝜽 and the empirical quantiles of order
hypothesis of convergence given in Eq. (12) is satisfied, therefore DM 𝐵 𝑣=1 𝑃 𝐵,𝑣
∑ 𝑞
𝑞, 0 < 𝑞 < 1, 𝑃̂ ̂ ∗ . These quantiles give an image of the scatter of
can be applied. Thus when 𝑁 = 𝑚 𝑗=1 𝑛𝑗 is large enough, 𝑔(𝜽̂ 𝑀𝐿 ) follows
( ) ∗
𝜽𝑃 𝐵
the normal distribution:  𝑔(𝜽), 𝐽𝑔 (𝜽)𝐽 −1 (𝜽̂ 𝑀𝐿 )𝐽𝑔 (𝜽)𝑇 , where 𝐽𝑔 (𝜽) is 𝜽̂ 𝑃 𝐵 compared with 𝜽̂ 𝑀𝐿 and allow to derive order of magnitude
the Jacobian of 𝑔, i.e. the 𝑟×𝑝 = 1×3 matrix of partial derivatives of the of the statistical estimation uncertainty on 𝜽.
entries of 𝑔 with respect to the entries of 𝜽 and 𝐽 is defined in Eq. (13).
For this purpose, considering the Birnbaum-Saunders approximation We will consider 𝜽̂ 𝑃 𝐵 is the natural bootstrap estimator of 𝜽, with
∗ 1 ∑𝐵 ̂ ∗
given in (11), we need to invert the Birnbaum-Saunders CDF and 𝜽̂ 𝑃 𝐵 = 𝜽̂ 𝑃 𝐵 = 𝜽 .
𝐵 𝑣=1 𝑃 𝐵,𝑣

6
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

3.3.2. Non-parametric bootstrap Table 1

Statistical indicators used to assess bootstrap methods performance.
Non-parametric bootstrap techniques consist in resampling data like
Mathematical expression
parametric bootstrap method except that data are not assumed to
∗
belong to a given parametric model. Bias 𝜽̂ − 𝜽̂ 𝑀𝐿
̂𝜽∗ − 𝜽̂ 𝑀𝐿
Moving block bootstrap (MBB). When data are organized in time series, Relative bias
𝜽̂ 𝑀𝐿
resampling single observations as considered by Efron and Tibshirani ( )2
1 ∑𝐵 ∗ ∗
Variance 𝜽̂ 𝑣 − 𝜽̂
[48] can fail to produce valid approximations in presence of depen- 𝐵 𝑣=1
√
dence (in time or between components) [54]. Thus, to solve this variance
Relative error
issue, Moving Block Bootstrap (MBB) can be suitable to efficiently 𝜽̂
∗

resample time-dependent data such as erosive cavitation successive Dispersion (at 95%) 𝑃̂ ̂97.5%
∗ − 𝑃̂ ̂2.5%
∗
𝜽 𝜽
measurements. This resampling method relies on retaining neighbors Mean squared error bias2 + variance
observations to form blocks. This way, the dependence structure of data
at short distances is preserved, Kreiss and Lahiri [54].
Let us consider the 𝑗th degradation path of length 𝑛𝑗 of a gamma
be considered as the EB estimate of 𝜽, and the empirical quantiles of
process as a time series and that degradation increments 𝛿𝑖,𝑗 are ob-
order 𝑞, 0 < 𝑞 < 1, 𝑃̂ ̂𝑞∗ .
served between discrete times 𝑡𝑖,𝑗 and 𝑡𝑖+1,𝑗 , 0 ≤ 𝑖 ≤ 𝑛𝑗 − 1. We denote 𝜽𝐸𝐵
𝐘𝑗 = {𝑌0,𝑗 , … , 𝑌𝑛𝑗 −1,𝑗 }, with 𝑌𝑖,𝑗 = {𝑡𝑖,𝑗 , 𝑡𝑖+1,𝑗 , 𝛿𝑖,𝑗 }, 1 ≤ 𝑗 ≤ 𝑚. Let 𝑙 be an
3.3.3. Bootstrap uncertainty on hitting time 𝑇𝜌
integer with 1 ≤ 𝑙 ≤ 𝑛𝑗 . The overlapping blocks 𝐵0,𝑗 , … , 𝐵𝑝,𝑗 , 𝑝 = 𝑛𝑗 − 𝑙 ∗
From the 𝐵 obtained 𝜽̂ 𝑣 , 1 ≤ 𝑣 ≤ 𝐵, from previous bootstrap
of length 𝑙 contained in 𝐘𝑗 are defined as follow:
techniques, the percentiles (of order 𝑞) of first hitting time distribution,
∗
𝐵0,𝑗 = (𝑌0,𝑗 , 𝑌1,𝑗 , … , 𝑌𝑙−1,𝑗 ), noted F−1
𝐵𝑆
(𝑞|𝜽̂ 𝑣 ), can be easily derived using the Birnbaum-Saunders
approximation (see Section 3.2).
𝐵1,𝑗 = (𝑌1,𝑗 , … , 𝑌𝑙−1,𝑗 , 𝑌𝑙,𝑗 ),
... ... 3.4. Statistical indicators for performance assessment
𝐵𝑝,𝑗 = (𝑌𝑛𝑗 −𝑙,𝑗 , … , 𝑌𝑛𝑗 −1,𝑗 ).
In Table 1, we introduce some classic statistical indicators to assess
For state of simplicity, consider that 𝑛𝑗 is a multiple of 𝑙. Let 𝑏𝑗 =
bootstrap techniques performance and compare them with MLE + DM.
𝑛𝑗 ∕𝑙. To generate one MBB sample, 𝑏𝑗 blocks are taken at random
If we consider the case of model parameters’ estimation, bias measures
with replacement from the collection {𝐵0,𝑗 , … , 𝐵𝑝,𝑗 }. The bootstrapped ̂ Sim-
∗ ,…,𝑌∗ how well the true unknown parameter 𝜽 is approximated by 𝜽.
dataset 𝐘𝑗 = {𝑌0,𝑗 } contains 𝑏𝑗 ⋅ 𝑙 = 𝑛𝑗 elements since each
𝑛𝑗 −1,𝑗 ilarly, variance and dispersion give an indication on the estimation
resampled block has 𝑙 elements and 𝑏𝑗 blocks are drawn. For time- accuracy of 𝜽 by 𝜽.̂ Finally, the Mean Squared Error (MSE) combines
1∕𝑘
dependent data, 𝑙 should satisfy the condition: 𝑙 = 𝐶𝑛𝑗 , with 𝑘 = 3 or both indicators and allows to compare techniques between each other.
4, where 𝐶 ∈ R is a constant [54]. In the general case, where 𝑙 does not ∗ ∗
𝜽̂ denotes the mean of 𝜽̂ , 1 ≤ 𝑣 ≤ 𝐵, as introduced for each bootstrap
𝑣
divide 𝑛𝑗 , 𝑏𝑗 = 𝑏0 blocks may be resampled with 𝑏0 = min{𝑘 ≥ 1 ∶ 𝑘𝑙 ≥
technique in the previous sections.
𝑛𝑗 } and only the first 𝑛𝑗 resampled data-values are retained to define
the bootstrap observations 𝐘∗𝑗 . The following algorithm illustrates the
4. Case study
MBB procedure:
4.1. Experimental dataset
1. Set the block length 𝑙.
2. Draw with replacement a block of 𝑙 consecutive observations
In this study, we aim to model erosive cavitation using NHGP
from 𝐘𝑗 .
and evaluate confidence intervals of model parameters with limited
3. Repeat 𝑏𝑗 times step 2. until the bootstrap sample 𝐘∗𝑗 has 𝑏𝑗 ⋅ 𝑙 ≥
dataset. In [55], a sensitivity analysis was conducted on gamma process
𝑛𝑗 elements.
parameters and size of dataset using Monte Carlo simulations. The latter
4. Retain only the first 𝑛𝑗 resampled values. study shows the performances of maximum likelihood estimation vary
5. Repeat steps 2–4 for each degradation path 𝑗, 1 ≤ 𝑗 ≤ 𝑚, to get depending on the degradation acceleration (parameter 𝑏 of the shape
a set of 𝑚 paths 𝐘∗ = {𝐘∗1 , … , 𝐘∗𝑚 }. function) and the number of paths 𝑚. In particular, authors show a
6. Repeat steps 5. 𝐵 times to generate 𝐵 sets of bootstrap samples dataset with more degradation paths should be preferred for a same
𝐘∗𝑣 , 1 ≤ 𝑣 ≤ 𝐵. total number of observations 𝑁.
∗ Although we focused in the present paper on modeling cavitation
7. Compute 𝜽̂ 𝑀𝐵𝐵,𝑣 using the ML method for each bootstrap sam-
∗
ples 𝐘𝑣 . observed on hydraulic runner, the same methodology can be applied to
other cumulative and monotonic deterioration measured on other types
∗
Similarly to PB method, from the 𝐵 obtained 𝜽̂ 𝑀𝐵𝐵,𝑣 , 1 ≤ 𝑣 ≤ 𝐵, we of systems. Moreover, resulting hitting time distribution variability will
∗ 1 ∑ 𝐵 ∗ be computed according to the results from the previous confidence
can consider 𝜽̂ 𝑀𝐵𝐵 = 𝜽̂ 𝑀𝐵𝐵 = 𝜽̂ the natural MBB estimate
𝐵 𝑣=1 𝑀𝐵𝐵,𝑣 intervals (using DM and bootstrap techniques). We used a dataset 𝐃𝑜𝑏𝑠
of 𝜽, and the empirical quantiles of order 𝑞, 0 < 𝑞 < 1, 𝑃̂ ̂𝑞∗ . In the consisting in 𝑚 = 3 degradation paths presented in Fig. 3 (right)
𝜽𝑀𝐵𝐵
case study presented in the next section, the block length 𝑙 was chosen obtained from a vibratory cavitation test (2 replications) performed
equal to 3 to make a compromise between time dependence structure in laboratory conforming to standard given in [38]. During test 𝑗, the
and the relative short path length 𝑛𝑗 . degradation level, here the specimen mass, is measured at instants 𝑡𝑖,𝑗 ,
1 ≤ 𝑖 ≤ 𝑛𝑗 , to determine the cumulative mass loss 𝑥𝑖,𝑗 as depicted
Efron bootstrap (EB). The bootstrap method originally introduced by
by degradation path 𝑗 in Fig. 1. Meanwhile, for actual Francis turbine
Efron [47] can be seen as a special case of MBB. Indeed if we set
prototypes, the incubation and acceleration regions of the erosive
𝑙 = 1, the MBB is equivalent to the ordinary non parametric bootstrap
cavitation pattern (shown in Fig. 3 (left)) are not observed. For this
method described by Efron and Tibshirani [48] for i.i.d. data. Similarly
∗ reason, we decided to only focus on modeling the degradation on the
to previous methods, from the 𝐵 obtained 𝜽̂ 𝐸𝐵,𝑣 , 1 ≤ 𝑣 ≤ 𝐵, we can
maximum rate region (ie: from time = 𝑡𝐴,𝑗 to 𝑡𝑛𝑗 ). First, degradation
∗ 1 ∑𝐵 ̂ ∗
introduce the empirical mean 𝜽̂ 𝐸𝐵 = 𝜽̂ 𝐸𝐵 = 𝜽 which will data in the region of interest have been transformed to satisfy the first
𝐵 𝑣=1 𝐸𝐵,𝑣

7
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

Table 2 Table 3
Estimations of 𝜽 for ML method and 3 different bootstrap techniques - Approximate Statistical indicators from different methods’ results.
asymptotic confidence interval (at 95%) of 𝜽. ML Bootstrap techniques
𝑢̂ 𝑐̂ 𝑏̂
PB EB MBB
ML 32.501 1.2722 1.1348
𝑢̂ 2.653 2.470 3.919
Approximate asymptotic CI (at 95%) [16.32, 48.69] [1.101, 1.939] [1.108, 1.162]
Bias 𝑐̂ 0 (supposed) 0.6794 0.0850 0.0827
PB 35.154 1.9516 1.0888 𝑏̂ −0.0461 0.0013 0.0073
EB 34.971 1.3572 1.1361
𝑢̂ 0.0816 0.0760 0.1206
MBB 36.420 1.3549 1.1421
Relative bias 𝑐̂ 0 (supposed) 0.5340 0.0668 0.0650
𝑏̂ −0.0406 0.0011 0.0064
𝑢̂ 68.20 87.56 48.59 58.11
Variance 𝑐̂ 0.1158 0.2926 0.0787 0.0865
gamma process property, 𝑋(0) = 0, with 𝑡′𝑖,𝑗 = 𝑡𝑖,𝑗 − 𝑡𝐴,𝑗 and 𝑥′𝑖,𝑗 =
𝑏̂ 0.0002 0.0002 0.0002 0.0004
𝑥𝑖,𝑗 − 𝑥𝐴,𝑗 , 𝐴 ≤ 𝑖 ≤ 𝑛𝑗 . Note that 𝑡𝐴,𝑗 (considered to be deterministic) is
𝑢̂ 0.2541 0.2879 0.2145 0.2346
determined by the intersection between the 𝑥-axis and the slope of the Relative error 𝑐̂ 0.2675 0.4252 0.2206 0.2312
degradation path in the maximum stage region. Further explanations on 𝑏̂ 0.0120 0.0109 0.0108 0.0182
the standard test method can be found, as detailed in [38]. Because we 𝑢̂ 32.64 36.77 25.84 29.31
notice in Fig. 3 (right) an increasing trend in degradation over time, we Dispersion at 95% 𝑐̂ 0.838 2.099 1.040 1.133
chose to model it with a NHGP (a performance comparison, conducted 𝑏̂ 0.0540 0.0486 0.0485 0.0814
previously in [56] on the same dataset, between homogeneous and 𝑢̂ 68.195 94.69 56.74 73.89
non-homogeneous gamma process has shown better results for NHGP). MSE 𝑐̂ 0.1158 0.7544 0.08961 0.09553
Thus, 𝜽̂ 𝑀𝐿 = (𝑢̂ 𝑀𝐿 , 𝑐̂𝑀𝐿 , 𝑏̂ 𝑀𝐿 ) was computed from the 𝑚 = 3 paths of 𝑏̂ 0.0002 0.0023 0.0002 0.0003
∑
length 𝑛𝑗 ∈ {10, 10, 14}, 1 ≤ 𝑗 ≤ 3, for a total of 𝑁 = 𝑚=3 𝑗=1 𝑛𝑗 = 34
observations, using ML equations (see Section 2.4) and optimization Table 4
tools. Because of numerical difficulties induced by the gamma function, Estimation of hitting time percentiles: 𝑃𝑇2.5% , 𝑃𝑇50% , 𝑃𝑇97.5% and their associated estimation
𝜌 𝜌 𝜌

Eqs. (5) to (7) have been used instead of maximizing directly expression uncertainty for the 3 bootstrap methods and ML + DM.

(4). Quartiles of 𝑇𝜌 ML+DM PB EB MBB

Median 976 950 976 976
4.2. Calibration of bootstrap sample number 𝑃𝑇2.5% Dispersion at 95% 31 40 40 54
𝜌
Relative error 0.032 0.042 0.041 0.055

In order to determine a relevant value of 𝐵, the number of bootstrap Median 1006 979 1005 1006
𝑃𝑇50% Dispersion at 95% 31 40 40 54
datasets to be drawn (see Section 3.3.2), a convergence analysis was 𝜌
Relative error 0.031 0.041 0.040 0.054
performed. For this purpose, we considered the previous dataset 𝐃𝑜𝑏𝑠
Median 1037 1011 1036 1036
and the three bootstrap methods were applied for different values 𝑃𝑇97.5% Dispersion at 95% 34 44 41 56
of 𝐵, 10 ≤ 𝐵 ≤ 104 . Fig. 4 shows 𝜽̂ 𝐸𝐵 and percentiles, 𝑃̂ ̂2.5%
∗ and 𝜌
Relative error 0.033 0.044 0.040 0.054
𝜽𝐸𝐵
𝑃̂ ̂97.5%
∗ , respectively represented with upper and lower bars as defined in
𝜽𝐸𝐵
Section 3.3. We observe the asymptotic convergence is reached around
𝐵 = 103 draws. We retained 𝐵 = 104 draws for the three bootstrap close to zero. One of the main hypothesis of approximate asymptotic
techniques to ensure a good convergence. confidence intervals from MLE is highlighted in Table 3: 𝜽̂ 𝑀𝐿 is consid-
ered centered on true unknown parameter 𝜽 and hence shows no bias
4.3. Estimation uncertainties for model parameters (asymptotically, which may be a wrong hypothesis for small dataset).
Nonetheless, considering ML has no bias, the MSE is the smallest for EB
The three different bootstrap methods presented in Section 3.3 method on both 𝑢, 𝑏 and 𝑐 parameters. We observe in Fig. 5 boxplots
have been applied on the experimental dataset introduced in Sec- are asymmetrical while approximate asymptotic CI are centered around
∗
tion 4.1. Approximate asymptotic confidence intervals introduced in 𝜽̂ 𝑀𝐿 which is not confirmed by bootstrap distributions of 𝜽̂ . Tables 2,
Section 3.3.1 were derived to get the approximate symmetrical two- 3 and Fig. 5 show a small relative error on parameter 𝑏 estimations
sided 95% confidence intervals on 𝜽. (around 1%) for non-parametric methods and exhibit value of power
Table 3 presents statistical indicators from bootstrap techniques and 𝑏 = 1 is not included in confidence interval at 95% level, thus justifying
ML results. For ML method, results are derived from formulas (12) the use of NHGP to model our cavitation erosion dataset.
and (16). We observed in this study, that the best methods (EB and
MBB), with the relative error reaching around 25% for parameters 𝑢
and 𝑐, which can be considered as quite high. Meanwhile, EB seems to 4.4. Influence of model parameter uncertainties on hitting-time
produce better estimation of 𝜽 than MBB in terms of relative estimation
uncertainties. For better illustration and comparison, Fig. 5 shows
Fig. 6 illustrates the concept of propagation of uncertainties. Indeed,
boxplots of 𝜽̂∗𝑣 , 1 ≤ 𝑣 ≤ 𝐵, for the three bootstrap techniques. Lower
influence of 𝜽 estimation uncertainty on hitting time distributions F𝑇𝜌
and upper whiskers denote 𝑃̂ ̂2.5%∗ and 𝑃̂ ̂97.5%
∗ respectively; box shows
𝜽 𝜽 can be derived using the Birnbaum-Saunders approximation from both
first 𝑃 ∗ and third 𝑃 ∗ quartiles and median 𝑃̂ 50%
̂ 25% ̂ 75%
∗ . Left and right delta method (see Section 3.2) and bootstrap techniques results. To
𝜽̂ 𝜽̂ 𝜽̂
dotted lines respectively represent the lower and upper bounds of the study how uncertainties impact the whole distribution of 𝑇𝜌 (threshold
asymptotic two-sided 95% confidence intervals in Table 2, obtained by 𝜌 = 100 mg), we compute in Fig. 6 for each percentile order 𝑞 =
the asymptotic normality of MLE. The centered dotted line represents {2.5%, 5%, … , 97.5%}: F−1 𝐵𝑆
(𝑞|𝜽̂ 𝑀𝐿 ) and its variability using DM (see
∗
𝜽̂ 𝑀𝐿 = (𝑢̂ 𝑀𝐿 , 𝑐̂𝑀𝐿 , 𝑏̂ 𝑀𝐿 ). In both Table 3 and Fig. 5 we observe the −1 ̂
Section 3.2) and F (𝑞|𝜽 ), with 1 ≤ 𝑣 ≤ 𝐵, for the three bootstrap
𝐵𝑆 𝑣
largest (positive, leading to a possible over-estimation) bias (but also techniques (see Section 3.3). In Fig. 6, dotted lines and shaded areas
∗
relative bias) for PB method (for parameters 𝑐 and 𝑏, around 10% and represent respectively medians and 𝑃̂ 75%𝑞 − 𝑃̂ 25%
𝑞 of F−1 (𝑞|𝜽̂ 𝑣 ). We
𝑃𝑇 𝐵𝑆 𝑃𝑇
50% respectively) and the largest scatter on 𝑢 and 𝑐, while EB and 𝜌 𝜌

MBB bias is around 10%. In terms of bias and scatter, EB and MBB observe clearly that the dispersion seems to remain constant for the
show similar results except for EB estimation of 𝑏 which exhibits a bias different percentile orders 𝑞 for all bootstrap methods. As observed in

8
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

Fig. 4. 𝜽̂ 𝐸𝐵 , 𝑃̂ ̂2.5%
∗ and 𝑃̂ ̂97.5%
∗ - Convergence analysis on bootstrap iteration number 𝐵 for EB method.
𝜽𝐸𝐵 𝜽𝐸𝐵

Fig. 5. Comparison of bootstrap and asymptotic CI.

Fig. 5 with boxplots, the shaded areas in Fig. 6 are slightly asymmetri-
cal compared to median (more significantly for MBB technique): most
values are concentrated on lower half of distribution. ML + DM, EB and
MBB produce very close results on both percentiles of lower and upper
orders but also on percentile estimation of central order. For ML + DM
and EB methods, we observe also close results in terms of dispersion of
percentile estimations. Only MBB produces larger dispersion on hitting
time percentile distributions. If we look closely at the relative error
of 𝑃𝑇𝑞 in Table 4, we can stress that 𝑇𝜌 seems to be better estimated
𝜌
than 𝜽. This result could be explained by an offset effect of estimation
uncertainty of parameters 𝑢, 𝑐 and 𝑏 on the hitting time distribution.

5. Discussion

In the light of results presented in Section 4.3, the use of a NHGP

can be justified considering CI of parameter 𝑏 does not include value 1,
i.e. the homogeneous case. Non-parametric bootstrap techniques seem
to provide results close to ML method, while parametric bootstrap
presents both the largest bias and scatter on parameter estimation.
Even if this may seem counter-intuitive, we noticed PB tends to fail
Fig. 6. Hitting time distributions of 𝑇𝜌 (for 𝜌 = 100 mg) for 3 bootstrap methods and
producing ‘‘good’’ simulated paths when their length 𝑛𝑗 is small (lower ML + DM.
than 103 increments). This trend is particularly true for NHGP, with a
MLE convergence performance worse than HGP for short path length.
MBB and EB give a good estimation of the degradation acceleration methods have more difficulties when it comes to estimate these param-
depicted by power parameter 𝑏, with relative bias and dispersion lower eters. Among bootstrap techniques, EB seems to produce the closest
than 1% for both methods. The same indicator reaches around 10% estimations relatively to 𝜽̂ 𝑀𝐿 . This result can be surprising considering
for parameters 𝑢 and 𝑐, which denotes that non parametric bootstrap

9
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

MBB is supposed to better capture temporal dependence of time series. influence on degradation phenomenon of covariates, such as material
This may be explained by the relatively small dataset (short path properties of studied systems or even operating conditions, to develop
length) used to resample; the new resampled dataset does not provide specific sub-models. These improvements may lead to more realistic
much more information compared to the original one. Note that this degradation prediction and maintenance planning.
remark applies also for other bootstrap techniques.
According to Table 3 and Fig. 6, we observe that the estimation
CRediT authorship contribution statement
uncertainty on parameter 𝑏 has the greatest influence on hitting time
distribution. The larger the scatter on parameter 𝑏 is, the larger 𝑇𝜌 per-
Q. Chatenet: Developed the theory, Performed the computations,
centile estimate scatter at a given order will be. A further investigation
Investigate delta method, Supervised the findings of this work, Dis-
should be carried out to quantify the influence of each parameter on
cussed the results, Contributed to the final manuscript. E. Remy:
hitting time 𝑇𝜌 using a sensitivity analysis. Combination of ML and DM
Verified the analytical methods, Discussed the results, Contributed to
seems to give results close to EB, for a computing lower cost: around
the final manuscript. M. Gagnon: Provided industrial dataset, Dis-
0.0029 s1 for ML + DM and around 2035 s for bootstrap techniques
cussed the results, Contributed to the final manuscript. M. Fouladirad:
(to generate 𝐵 = 104 datasets and estimate the associated parameters).
Verified the analytical methods, Discussed the results, Contributed to
Meanwhile, mathematical developments require much more effort for
the final manuscript. A.S. Tahan: Discussed the results, Contributed to
the DM to obtain uncertainty on hitting time distribution. Moreover,
the final manuscript.
bootstrap techniques give the entire distribution of parameters’ es-
timators instead of only brackets centered on 𝜽̂ 𝑀𝐿 for approximate
asymptotic confidence intervals for ML method. This feature allows, Declaration of competing interest
in the case of bootstrap techniques, to estimate possible bias between
𝜽̂ 𝑀𝐿 and 𝜽. For these reasons and because ML relies on asymptotic The authors declare that they have no known competing finan-
hypothesis (not necessarily satisfied), EB shows better performances cial interests or personal relationships that could have appeared to
and tends to be more robust than the combination of ML and DM. influence the work reported in this paper.

6. Conclusion and perspectives

Acknowledgment
In this paper we modeled erosive cavitation from laboratory exper-
iments with a stochastic process taking into account physical random The work was supported by MITACS.
dynamics. Maximum likelihood method was used to estimate the un-
known parameters of the NHGP. Moreover, several methods were used Appendix
to assess model parameters estimation uncertainties and their influence
on the distribution of the hitting time for a given degradation level
For a gamma process with shape function 𝑣(𝑡) = 𝑐𝑡𝑏 , first and
threshold. To our knowledge, we are the first to apply delta method
second partial derivatives of log-likelihood function 𝓁(𝜽), respectively
in the framework of gamma processes. It is found on our case that
to parameter 𝜽, are given in the following expressions:
the combination of ML + DM, and Efron bootstrap produce similar
results. Meanwhile, bootstrap methods exhibit quicker and simpler 𝜕𝓁(𝜽) ∑ ∑𝑛𝑗 𝑐 [ 𝑏 ]
𝜕𝑢
= 𝑚 𝑗=1 𝑖=1 𝑢 𝑖,𝑗
𝑡 − 𝑡𝑏𝑖−1,𝑗 − 𝛿𝑖,𝑗
implementation than delta method. Similarly, MBB shows a reasonable
𝜕𝓁(𝜽) ∑ ∑ 𝑛𝑗 [ 𝑏 ][ ( )]
bias but more dispersion. On the contrary, PB fails to give ‘‘accurate’’
𝜕𝑐
= 𝑚𝑗=1 𝑡 − 𝑡𝑏𝑖−1,𝑗 log(𝛿𝑖,𝑗 ) + log(𝑢) − 𝜓 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ]
𝑖=1 𝑖,𝑗
estimates of parameters. However, reader needs to keep in mind the
𝜕𝓁(𝜽) ∑𝑚 ∑𝑛𝑗 [ 𝑏 𝑏
]
reference dataset show a particularly small size and general conclusions 𝜕𝑏
= 𝑗=1 𝑖=1
𝑐 𝑡 𝑖,𝑗 log(𝑡 𝑖,𝑗 ) − 𝑡 𝑖−1,𝑗
log(𝑡 𝑖−1,𝑗 )
should not be drawn from this specific case. In particular, data were [ ( ) ]
× log(𝛿𝑖,𝑗 ) − 𝜓 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] + log(𝑢)
collected during experiments under controlled environment showing
small discrepancy. Meanwhile, random effect models can be useful to (A.1)
consider unexplained heterogeneous degradation rates within a same
product population or a same test procedure [23,57]. For a further ∑𝑚 ∑𝑛𝑗 [ ]
𝜕 2 𝓁(𝜽) 𝑐
𝜕𝑢2
=− 𝑗=1 𝑖=1 𝑢2
𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗
study on the performances of ML parameter estimation, we refer the
reader to [55]. 𝜕 2 𝓁(𝜽) ∑ ∑𝑛𝑗 [ 𝑏 ]2 ( )
𝜕𝑐 2
=− 𝑚 𝑗=1 𝑡 − 𝑡𝑏𝑖−1,𝑗 𝜓 ′ 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ]
𝑖=1 𝑖,𝑗
Because of the small dataset and the potential availability of other
𝜕 2 𝓁(𝜽) ∑𝑚 ∑𝑛𝑗 {[ 𝑏 2 − 𝑡𝑏 2
]
sources of information about cavitation erosion, such as expert knowl- 2 = 𝑗=1 𝑖=1
𝑐 𝑡 𝑖,𝑗 log(𝑡 𝑖,𝑗 ) 𝑖−1,𝑗
log(𝑡 𝑖−1,𝑗 )
𝜕𝑏
edge, monitoring data or results from numerical simulations, a Bayesian [ ( ) ]
framework should be well suited to combine these multiple sources of × log(𝑡𝑖,𝑗 ) − 𝜓 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] + log(𝑢)
[ ]2 [ ( )]
information and decrease uncertainty on parameters estimations [58, −𝑐 𝑡𝑏𝑖,𝑗 log(𝑡𝑖,𝑗 ) − 𝑡𝑏𝑖−1,𝑗 log(𝑡𝑖−1,𝑗 ) 𝜓 ′ 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] }
59]. Further investigation could be carried out using other stochastic
𝜕 2 𝓁(𝜽) ∑𝑚 ∑𝑛𝑗 1 [ 𝑏 𝑏
]
processes or degradation models such as degradation paths [11], cu- = 𝑗=1 𝑖=1
𝑡𝑖,𝑗 − 𝑡𝑖−1,𝑗
𝜕𝑢𝜕𝑐 𝑐
mulative shock models [30] or inverse gaussian process [23]. Likewise, 𝜕 2 𝓁(𝜽) ∑𝑚 ∑𝑛𝑗 𝑐 [ 𝑏 𝑏
]
(A.2)
other bootstrap techniques, such as Markov bootstrap [60], could also 𝜕𝑢𝜕𝑏
= 𝑗=1 𝑖=1 𝑣 𝑖,𝑗
𝑡 log(𝑡 𝑖,𝑗 ) − 𝑡𝑖−1,𝑗 log(𝑡𝑖−1,𝑗 )
be investigated. 𝜕 2 𝓁(𝜽) ∑ ∑𝑛𝑗 1 [ 𝑏 ]
𝜕𝑐𝜕𝑢
= 𝑚 𝑗=1 𝑖=1 𝑢 𝑖,𝑗
𝑡 − 𝑡𝑏𝑖−1,𝑗
This paper can be seen as the first step in development of condition
𝜕 2 𝓁(𝜽) ∑ ∑𝑛𝑗 [ 𝑏 ]
based maintenance policies for high cost systems as presented in [61]
𝜕𝑐𝜕𝑏
= 𝑚 𝑗=1 𝑡 log(𝑡𝑖,𝑗 ) − 𝑡𝑏𝑖−1,𝑗 log(𝑡𝑖−1,𝑗 )
𝑖=1 𝑖,𝑗
and Lei et al. [62]. Although this study focuses on cavitation ero- ( )
sion, a similar methodology could be applied to other cumulative and × {log(𝛿𝑖,𝑗 ) − 𝜓 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] + log(𝑢)
( )
monotonic degradation observed on different systems. Next steps would −𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ]𝜓 ′ 𝑐[𝑡𝑏𝑖,𝑗 − 𝑡𝑏𝑖−1,𝑗 ] }
imply using real monitoring cavitation measurements and study the
𝜕 2 𝓁(𝜽) 𝜕 2 𝓁(𝜽)
𝜕𝑏𝜕𝑢
= 𝜕𝑢𝜕𝑏
𝜕 2 𝓁(𝜽) 𝜕 2 𝓁(𝜽)
1
Performed on a quad core Intel i7 2.8 GHz processor with 16 GB RAM. 𝜕𝑏𝜕𝑐
= 𝜕𝑐𝜕𝑏

10
Q. Chatenet et al. Reliability Engineering and System Safety 213 (2021) 107671

References [35] Nystad BH, Gola G, Hulsund JE, Roverso D. Technical condition assessment and
remaining useful life estimation of choke valves subject to erosion. In: Annual
[1] Nguyen KT, Fouladirad M, Grall A. Model selection for degradation modeling and Conference of the Prognostics and Health Management Society. 2010, p. 11–3.
prognosis with health monitoring data. Reliab Eng Syst Saf 2018;169:105–16. [36] Arfken G. The incomplete gamma function and related functions. Math Methods
[2] Verma AK, Ajit S, Karanki DR. Reliability and safety engineering, Vol. 43. Phys 1985;565–72.
Springer; 2010. [37] Park C, Padgett W. Accelerated degradation models for failure based on
[3] Kahle W, Mercier S, Paroissin C. Degradation processes in reliability. John Wiley geometric brownian motion and gamma processes. Lifetime Data Anal
& Sons; 2016. 2005;11(4):511–27.
[4] Escaler X, Egusquiza E, Farhat M, Avellan F, Coussirat M. Detection of cavitation [38] Standard Test Method for Cavitation Erosion Using Vibratory Apparatus. West
in hydraulic turbines. Mech Syst Signal Process 2006;20(4):983–1007. Conshohocken, PA: ASTM International; 2016, www.astm.org.
[5] Brennen CE. Cavitation and bubble dynamics. Cambridge University Press; 2013. [39] Paroissin C, Salami A. Failure time of non homogeneous gamma process. Comm
[6] Celebioglu K, Altintas B, Aradag S, Tascioglu Y. Numerical research of cavitation Statist Theory Methods 2014;43(15):3148–61.
on Francis turbine runners. Int J Hydrogen Energy 2017;42(28):17771–81. [40] Balakrishnan N, Kundu D. Birnbaum-saunders distribution: A review of models,
[7] Krumenacker L, Fortes-Patella R, Archer A. Numerical estimation of cavitation analysis and applications. 2018, arXiv preprint arXiv:1805.06730.
intensity. In: IOP Conference Series: Earth and Environmental Science, Vol. 22. [41] Birolini A. Reliability engineering, Vol. 5. Springer; 2007.
(5). IOP Publishing; 2014. [42] Hoadley B. Asymptotic properties of maximum likelihood estimators for the
[8] Gohil P, Saini R. Numerical study of cavitation in francis turbine of a small independent not identically distributed case. Ann Math Stat 1971;1977–91.
hydro power plant. J Appl Fluid Mech 2016;9(1). [43] Rao CR. Linear statistical inference and its applications, Vol. 2. Wiley New York;
[9] Singpurwalla ND. Survival in dynamic environments. Stat Sci 1995;86–103. 1973.
[10] Elsayed EA. Reliability engineering. John Wiley & Sons; 2020. [44] Efron B, Hinkley DV. Assessing the accuracy of the maximum likeli-
[11] Xu Z, Hong Y, Jin R. Nonlinear general path models for degradation data with hood estimator: Observed versus expected Fisher information. Biometrika
dynamic covariates. Appl Stoch Models Bus Ind 2016;32(2):153–67. 1978;65(3):457–83.
[12] Van Noortwijk J. A survey of the application of gamma processes in maintenance. [45] Casella G, Berger RL. Statistical inference, Vol. 2. Duxbury Pacific Grove, CA;
Reliab Eng Syst Saf 2009;94(1):2–21. 2002.
[13] Alaswad S, Xiang Y. A review on condition-based maintenance optimiza- [46] Fouladirad M, Paroissin C, Grall A. Sensitivity of optimal replacement policies
tion models for stochastically deteriorating system. Reliab Eng Syst Saf to lifetime parameter estimates. European J Oper Res 2018;266(3):963–75.
2017;157:54–63. [47] Efron B. Bootstrap methods: another look at the jackknife. In: Breakthroughs in
[14] Grinstead CM, Snell JL. Introduction to probability. American Mathematical Soc.; Statistics. Springer; 1992, p. 569–93.
2012. [48] Efron B, Tibshirani RJ. An introduction to the bootstrap. CRC press; 1994.
[15] Karlin S. A first course in stochastic processes. Academic press; 2014. [49] Vittinghoff E, Glidden DV, Shiboski SC, McCulloch CE. Regression methods in
[16] Kadloor S, Adve RS, Eckford AW. Molecular communication using brownian biostatistics: linear, logistic, survival, and repeated measures models. Springer
motion with drift. IEEE Trans NanoBiosci 2012;11(2):89–99. Science & Business Media; 2011.
[17] Zhang Z, Si X, Hu C, Lei Y. Degradation data analysis and remaining useful life [50] Kaneishi T, Dohi T. Parametric bootstrapping for assessing software reliability
estimation: A review on Wiener-process-based methods. European J Oper Res measures. In: Dependable Computing (PRDC), 2011 IEEE 17th Pacific Rim
2018. International Symposium on. IEEE; 2011, p. 1–9.
[18] Çınlar E. Probability and stochastics, Vol. 261. Springer Science & Business [51] Stute W, Manteiga WG, Quindimil MP. Bootstrap based goodness-of-fit-tests.
Media; 2011. Metrika 1993;40(1):243–56.
[19] Wang X, Xu D. An inverse Gaussian process model for degradation data. [52] Genest C, Rémillard B, et al. Validity of the parametric bootstrap for goodness-
Technometrics 2010;52(2):188–97. of-fit testing in semiparametric models. In: Annales de L’Institut Henri Poincaré,
[20] Abdel-Hameed M. A gamma wear process. IEEE Tans Reliab 1975;24(2):152–3. ProbabilitéS et Statistiques, Vol. 44. (6):Institut Henri Poincaré; 2008, p.
[21] Guan Q, Tang Y, Xu A. Objective Bayesian analysis accelerated degradation test 1096–127.
based on Wiener process models. Appl Math Model 2016;40(4):2743–55. [53] Avramidis A, L’ecuyer P, Tremblay P-A. Efficient simulation of gamma and
[22] Whitmore G. Estimating degradation by a Wiener diffusion process subject to variance-gamma processes, In: Proceedings of the 2003 Winter Simulation
measurement error. Lifetime Data Anal 1995;1(3):307–19. Conference, Vol. 1, 2003: p. 319–326.
[23] Peng C-Y. Inverse Gaussian processes with random effects and explanatory [54] Kreiss J-P, Lahiri SN. Bootstrap methods for time series. In: Handbook of
variables for degradation data. Technometrics 2015;57(1):100–11. Statistics, Vol. 30. Elsevier; 2012, p. 3–26.
[24] Ling MH, Tsui KL, Balakrishnan N. Accelerated degradation analysis for [55] Chatenet Q, Fouladirad M, Remy E, Gagnon M, Ton-That L, Tahan A. On the
the quality of a system based on the gamma process. IEEE Trans Reliab performance of the maximum likelihood estimation method for gamma process.
2015;64(1):463–72. In: European Safety and Reliability Conference, Vol. 25. 2019.
[25] Ye Z-S, Chen N. The inverse Gaussian process as a degradation model. [56] Chatenet Q, Gagnon M, Tôn-Thât L, Remy E, Fouladirad M, Tahan A. Stochas-
Technometrics 2014;56(3):302–11. tic modeling of cavitation erosion in francis runner. Int J Fluid Mach Syst
[26] Deng Y, Barros A, Grall A. Degradation modeling based on a time-dependent 2020;13(2):400–8.
ornstein-uhlenbeck process and residual useful lifetime estimation. IEEE Trans [57] Zhai Q, Chen P, Hong L, Shen L. A random-effects Wiener degradation model
Reliab 2016;65(1):126–40. based on accelerated failure time. Reliab Eng Syst Saf 2018;180:94–103.
[27] Castanier B, Grall A, Bérenguer C. A condition-based maintenance policy with [58] Bernardo JM, Smith AF. Bayesian theory, vol. 405. John Wiley & Sons; 2009.
non-periodic inspections for a two-unit series system. Reliab Eng Syst Saf [59] Bousquet N, Fouladirad M, Grall A, Paroissin C. Bayesian gamma processes for
2005;87(1):109–20. optimizing condition-based maintenance under uncertainty. Appl Stoch Models
[28] de Jonge B. Maintenance Optimization based on Mathematical Modeling. Bus Ind 2015;31(3):360–79.
University of Groningen; 2017. [60] Härdle W, Horowitz J, Kreiss J-P. Bootstrap methods for time series. Internat
[29] Deloux E, Fouladirad M, Bérenguer C. Health-and-usage-based maintenance Statist Rev 2003;71(2):435–59.
policies for a partially observable deteriorating system. Proc Inst Mech Eng Part [61] Zhu Q, Peng H, van Houtum G-J. A condition-based maintenance policy for
O 2016;230(1):120–9. multi-component systems with a high maintenance setup cost. Or Spectrum
[30] Zhu W, Fouladirad M, Bérenguer C. Condition-based maintenance policies for a 2015;37(4):1007–35.
combined wear and shock deterioration model with covariates. Comput Ind Eng [62] Lei Y, Li N, Guo L, Li N, Yan T, Lin J. Machinery health prognostics: A
2015;85:268–83. systematic review from data acquisition to RUL prediction. Mech Syst Signal
[31] Wang X. A pseudo-likelihood estimation method for nonhomogeneous gamma Process 2018;104:799–834.
process model with random effects. Statist Sinica 2008;18(3):1153–63.
[32] Le Son K, Fouladirad M, Barros A. Remaining useful lifetime estimation and
noisy gamma deterioration process. Reliab Eng Syst Saf 2016;149:76–87.
[33] Lu D, Pandey MD, Xie W-C. An efficient method for the estimation of parameters
of stochastic gamma process from noisy degradation measurements. Proc Inst
Mech Eng Part O 2013;227(4):425–33.
[34] Cholette ME, Yu H, Borghesani P, Ma L, Kent G. Degradation modeling and
condition-based maintenance of boiler heat exchangers using gamma processes.
Reliab Eng Syst Saf 2019;183:184–96.

11