Léa Brenière
email: lea.breniere@univ-grenoble-alpes.fr

Laurent Doyen
email: laurent.doyen@univ-grenoble-alpes.fr

Christophe Bérenguer
email: christophe.berenguer@grenoble-inp.fr

Optimization of Preventive Replacements Dates and Covariate Inspections for Repairable Systems in Varying Environments

Keywords: Reliability, Virtual age model, Imperfect maintenance, Dynamic covariate, Condition-based maintenance

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Technological or industrial equipment or structural components (e.g. dikes and dams) are subject to degradation because of intrinsic wear, usage imposed by operating conditions, or exposure to environmental factors. For such repairable industrial equipment, an important issue that was recently reviewed by de Jonge and Scarf (2020) is maintaining the system in working order conditions in accordance with safety, availability, and cost constraints. At present, as highlighted by [START_REF] Hong | Big data and reliability applications: The complexity dimension[END_REF], wide access to information and data on systems opens new research perspectives on this classical reliability and maintenance issue.

Two main approaches have been used in the reliability literature to address the above issues. One is based on recurrent event models, and is centred on probabilistic lifetime modelling of successive discrete event times and types (the successive times of failure, maintenance or replacement). Corresponding models can be uniquely defined through their conditional (failure) intensity, which generalises the failure time hazard rate as a conditional rate given the system history. The other considers stochastic processes (or even dynamic state-space models) that represent continuously time-varying degradation, in which the failures correspond to the first hitting time of a given threshold. Classical continuous degradation processes include Gamma processes, as reviewed by [START_REF] Van Noortwijk | A survey of the application of Gamma processes in maintenance[END_REF], and monotone Levy or Wiener processes, as in [START_REF] Zhang | Degradation data analysis and remaining useful life estimation: A review on Wiener-process-based methods[END_REF]. In particular, the discrete degradation models considered by [START_REF] Zhang | Continuous-observation partially observable semi-Markov decision processes for machine maintenance[END_REF] use Markov or semi-Markov chains and discretise the system health conditions into a finite number of states.

However, in many practical situations, the link between the degradation level and system failure is not direct, and the degradation does not lead automatically (with probability 1) to a failure when the degradation level reaches a given threshold. For example, this is the case when considering a crack length in a piece of structure, corrosion in a metallic part, or fatigue phenomena in mechanical parts. Hence, from a modelling perspective, such degradation/failure behaviour cannot be represented by a deterioration-threshold failure model, but rather, using a recurrent event model with a conditional failure intensity in which the degradation level is integrated as a covariate. This is the class of degradation/failure model that is used in this study. The above concept has already been used in the literature; see, for example, [START_REF] Liu | A dynamic prescriptive maintenance model considering system aging and degradation[END_REF], [START_REF] Zheng | Condition-based maintenance with dynamic thresholds for a system using the proportional hazards model[END_REF], [START_REF] Golmakani | Optimal replacement policy and inspection interval for condition-based maintenance[END_REF], [START_REF] Makis | Optimal replacement in the proportional hazards model[END_REF], [START_REF] Tian | Condition based maintenance optimization for multi-component systems using proportional hazards model[END_REF], [START_REF] Tian | Condition based maintenance optimization considering multiple objectives[END_REF], [START_REF] Vlok | Optimal component replacement decisions using vibration monitoring and the proportional-hazards model[END_REF], and [START_REF] Lam | A myopic policy for optimal inspection scheduling for condition based maintenance[END_REF]. Moreover, [START_REF] Liu | A dynamic prescriptive maintenance model considering system aging and degradation[END_REF] presented a real case application study for locomotive wheels of a Swedish company. According to the authors, the diameter of a wheel usually indicates the degradation level, where soft failure is defined when the diameter is reduced to a pre-specified threshold. However, during operation, the wheels are occasionally reprofiled owing to increased roughness, even if the wheel diameter remains above the tolerance level. Various factors affect the roughness of the wheels, including the weather conditions, surface smoothness of the track, and running speed, all of which may cause sudden failure even before the failure threshold is reached. Motivated by this limitation of degradationthreshold models, [START_REF] Zheng | Condition-based maintenance with dynamic thresholds for a system using the proportional hazards model[END_REF] proposed the integration of the degradation level effect into the failure time hazard rates.

Preventive maintenance actions and repairs may not all have the same effect on the system. Indeed, if the repair consists of a minimal intervention that is simply aimed at restarting the system without seeking to improve its condition, it will remain in the same state as directly before the failure. Such a repair action is termed as bad as old (ABAO); the repair is minimal and the failure time process is a non-homogeneous Poisson process (NHPP). Conversely, preventive maintenance action, such as a replacement that completely renews the system, is termed as good as new (AGAN), which is then perfect, and the failure time process is a renewal process. However, the reality may fall between these two situations, whereby the system directly after maintenance is better than old but worse than new and the maintenance is said to be imperfect. Many failure time process models exist to deal with imperfect maintenance, such as the Brown-Proschan model in [START_REF] Brown | Imperfect repair[END_REF], where the maintenance is ABAO or AGAN with a probability p or 1p. However, the present work relies on virtual age models, as exposed in [START_REF] Kijima | Some results for repairable systems with general repair[END_REF], and more precisely, on the parameterisation presented in [START_REF] Doyen | A generic framework for generalized virtual age models[END_REF], whereby the maintenance effect is characterised by an efficiency parameter ρ.

Many references can be found in the literature regarding the optimisation of preventive maintenance dates. These dates may follow a schedule that can be static and fixed beforehand or evolve dynamically throughout the life of the system. Among static schedules, periodic maintenance policies consider preventive replacements occurring at a fixed time interval, as presented in [START_REF] Barlow | Optimum preventive maintenance policies[END_REF]: preventive maintenance actions are periodic renewals of the system, and corrective maintenance actions are either minimal (ABAO) or perfect (AGAN). The aim is to determine the optimal replacement period τ P . Preventive maintenance actions are performed at times τ P , 2τ P , . . .. Many studies have assumed, as in [START_REF] Barlow | Optimum preventive maintenance policies[END_REF], that preventive maintenance actions are perfect. With minimal corrective maintenance actions, [START_REF] Gilardoni | Optimal maintenance time for repairable systems[END_REF] studied a periodic policy, including the estimation of the model parameters. In [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] and [START_REF] Zhang | A bivariate optimal repair-replacement model using geometric processes for a cold standby repairable system[END_REF], the corrective maintenance actions were imperfect. [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] proposed a dynamic policy that considered the failure process history to schedule the replacements. [START_REF] Zhang | A bivariate optimal repair-replacement model using geometric processes for a cold standby repairable system[END_REF] developed a model with two components that relieved one another following a failure, in which the repair durations were taken into account. The system is renewed either after a fixed number of failures or a certain working time. Conversely, several studies have investigated minimal corrective and imperfect preventive maintenance actions. For example, in [START_REF] Liao | Preventive maintenance scheduling for repairable system with deterioration[END_REF], preventive maintenance actions occur as soon as the system reliability reaches a threshold, and in [START_REF] Moghaddam | Preventive maintenance and replacement scheduling for repairable and maintainable systems using dynamic programming[END_REF], either a minimisation of the cost under a reliability constraint or a maximisation of the reliability under a budgetary constraint is applied. Finally, certain authors have considered imperfect corrective and preventive maintenance actions simultaneously. We can cite [START_REF] Nguyen | Analytical properties of an imperfect repair model and application in preventive maintenance scheduling[END_REF], who proposed three maintenance policies with the ARA ∞ -PLP model and relied on the properties of this model to solve the optimisation problem, or [START_REF] Safaei | A repair and replacement policy for repairable systems based on probability and mean of profits[END_REF], in which the dual goal was to find the optimal preventive maintenance period and to decide whether the system should be renewed or imperfectly repaired. Furthermore, we cite [START_REF] Love | A discrete semi-Markov decision model to determine the optimal repair/replacement policy under general repairs[END_REF], in which at each failure, the decision must be made to renew or imperfectly repair the system, but without considering preventive maintenance, and [START_REF] Zhou | Sequential imperfect preventive maintenance model with failure intensity reduction with an application to urban buses[END_REF], in which a periodic maintenance policy was outlined whereby the preventive maintenance actions have a variable effect, but without providing the optimal replacement period for a system.

The principle of condition-based maintenance is to integrate updated condition-monitoring information into the maintenance optimisation procedure. In the recurrent event modelling framework considered in this work, the proportional hazard model in [START_REF] Cox | Regression models and life tables (with discussion)[END_REF] in particular has been used in situations where the hazard rate of a system also depends on degradation processes such as erosion, corrosion, crack growth, and wear: these degradation processes (that can be modelled, for example, by Gamma processes or Markov chains) are treated as covariates and integrated into the hazard function through a link function. In almost all studies considering such condition-based maintenance policies, corrective maintenance actions are assumed to be perfect, as in [START_REF] Golmakani | Optimal replacement policy and inspection interval for condition-based maintenance[END_REF], [START_REF] Makis | Optimal replacement in the proportional hazards model[END_REF], [START_REF] Tian | Condition based maintenance optimization for multi-component systems using proportional hazards model[END_REF], [START_REF] Tian | Condition based maintenance optimization considering multiple objectives[END_REF], [START_REF] Vlok | Optimal component replacement decisions using vibration monitoring and the proportional-hazards model[END_REF] and [START_REF] Lam | A myopic policy for optimal inspection scheduling for condition based maintenance[END_REF]. In the final article, the covariate inspection cost was also considered. To the best of our knowledge, only [START_REF] Zheng | Condition-based maintenance with dynamic thresholds for a system using the proportional hazards model[END_REF] studied minimal corrective maintenance actions and time-dependent covariates (or degradation processes) in preventive replacement scheduling. However, no study has considered imperfect corrective maintenance actions in this context.

From a general perspective, covariates allow the heterogeneity among several independent replicas of identical systems to be accounted for (i.e., with the same failure model and the same parameters), whether or not this heterogeneity is observed. For example, [START_REF] Liu | Unobserved heterogeneity in stable imperfect repair models[END_REF] considered unobserved heterogeneity in virtual age models. In the case of observed heterogeneity, the covariates can either be fixed or time varying. Fixed static covariates can represent heterogeneity, such as the manufacturing origin or geographical location. Time-varying dynamic covariates can either integrate the effect of a degradation process (which can often be viewed as a composite degradation index) or directly variable operating conditions such as the temperature, pH, pressure, and workload. In this study, we focus on the dynamic covariates that represent degradation and we do not consider the case of variable operating conditions. In [START_REF] Brenière | Virtual age models with time-dependent covariates: A framework for simulation, parametric inference and quality of estimation[END_REF], a simulation method and an estimation procedure were established for virtual age models with observed fixed or timedependent covariates. Building on the same modelling framework, the present work investigates condition-based replacement policies; that is, how to use the information originating from different types of covariates to schedule preventive replacements. When developing a condition-based maintenance policy, the issue is not only to optimise the maintenance dates, but also to determine the optimal inspection time points at which covariates should be observed; for example, by considering the inspection cost in the formulation of the maintenance optimisation problem. Indeed, additional inspections would lead to better knowledge of the system, which should allow for more efficient scheduling of preventive maintenance actions. However, the gain in terms of maintenance costs may be exceeded by the cost of the inspections; thus, it is necessary to optimise the inspection and maintenance dates.

Thus, in this study, we develop a framework that aims to formulate and solve the optimisation of preventive replacement dates for systems that are influenced by covariates and undergo imperfect corrective repairs. We adopt a recurrent event modelling approach, and by extending the initial proposition of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], we propose a procedure to optimise preventive replacement dates by considering both the covariate information that is available online and imperfect corrective maintenance effects. The underlying principle of the proposed maintenance decision rule consists of determining the replacement time from a threshold criterion based on the conditional failure intensity of the maintained system. The case of static covariates is first investigated, and it is demonstrated how different levels of knowledge of the fixed covariates can be integrated into the maintenance decision and lead to different optimal replacement time points that are associated with different cost performances of the preventive replacement policy. A complete analysis of the design and performance of the different replacement policies is initially presented in the case of minimal repairs and then extended to the case of imperfect repairs. Subsequently, the proposed approach is generalised to condition-based maintenance, in which the dynamic covariates corresponding to the degradation process are integrated. We first examine the case of minimal repair and, in particular, the computation of the intensity threshold beyond which preventive maintenance must be triggered. Thereafter, we focus on optimising the covariate inspection period or date when the covariates are not completely known and must be inspected. Finally, a general case of imperfect repair is considered.

The remainder of this paper is organised as follows: Section 2 is dedicated to the background of this study. In particular, Sections 2.2 and 2.4 present the arithmetic reduction of age (ARA) model of [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF] and its generalisation to the case in which covariates are considered. The objective cost criteria are described in detail in Section 2.1. Section 2.3 summarises the work conducted by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] on the optimisation of preventive replacements with imperfect repairs, but without covariates, which is the starting point of our work. Section 3 analyses the case of fixed covariates, considering different levels of information on the covariates and minimal or imperfect repairs. Section 4 deals with the optimisation of preventive replacement dates with dynamic covariates and minimal repairs. It also presents a numerical procedure for the optimisation of the covariate inspection period and for the determination of the optimal date in the case of a single inspection. Section 5 extends the proposed approach to include imperfect repairs. Finally, Section 6 presents conclusions and provides insights for future work.

Background & problem statement

Preventive replacement policy and cost criterion

We consider a system undergoing imperfect corrective repairs after failure, which is eventually preventively replaced at the end of a maintenance cycle. Let us assume that the preventive replacements are perfect in the sense that they completely renew the system and incur a cost C P rev . These imperfect corrective repairs incur a cost C Corr . If C Corr represents the cost of a corrective maintenance action only; that is, with a lower efficiency than a perfect replacement, it is reasonably smaller than C P rev , which represents the total system renewal. If C Corr also considers all costs associated with a failure (such as unscheduled intervention and failure consequences), it can be higher than C P rev .

Let T 1 , . . . , T n denote the n first system failure times, with T 0 = 0, and {N t } t≥0 be the counting process of failures. Maintenance durations are not considered, and failure time points also correspond to corrective repair time points. Let {τ r } r≥0 be the successive replacement dates induced by preventive maintenance, with τ 0 = 0 (as the system is supposed to be new at the beginning of its life), {x r } r≥1 be the duration between two successive renewals; that is, for all r ≥ 1, x r = τ rτ r-1 , {R t } t≥0 is the counting process of replacements, and M r = C P rev + C Corr (N τr -N τ r-1 ), r ≥ 1 be the maintenance costs of each renewal cycle. Because the system is renewed by preventive replacements, the durations {x r } r≥1 and costs {M r } r≥1 are independent and have the same distributions. The total cost at time t ≥ 0 is equal to

C(t) = Rt r=1 M r = Rt r=1 (C P rev + C Corr (N τr -N τ r-1 )) + C Corr (N t -N τ R t ).
(1)

According to the renewal theorem (see, e.g. [START_REF] Tijms | A first course in stochastic models[END_REF]) with probability 1,

lim t→∞ C(t) t = EM EX , ( 2
)
where M is a random variable with the same distribution as {M r } r≥1 and X has the same distribution as {x r } r≥1 .

The preventive replacement date τ is the only decision variable of this preventive maintenance policy, which can be optimised by determining the value of τ that minimises the mean asymptotical cost per time unit:

G(τ ) = EM EX = E[C P rev + C Corr N τ ] Eτ = C P rev + C Corr EN τ Eτ .
(3)

Virtual age models

To model the corrective maintenance action efficiency; that is, imperfect repairs with intermediate effects between the ABAO and AGAN cases, we resort to virtual age models in this work.

The conditional failure intensity λ t of the counting process is defined as the conditional instantaneous failure rate:

λ t = lim ∆t→0 1 ∆t P(N t+∆t -N t -= 1 | F t -), (4)
where F t -and N t -are the process history and failure counting process directly before t, respectively. This representation characterizes the distribution of the next failure time given the process history. Let Λ t = t 0 λ s ds denote the cumulative intensity. In the virtual age model setting [START_REF] Kijima | Some results for repairable systems with general repair[END_REF], the system after a failure is in the same state as a new system that has survived until the virtual age without failure. This virtual age is denoted by V i (t), where t is the time and i is the number of failures. Following [START_REF] Doyen | A generic framework for generalized virtual age models[END_REF], a generic parameterisation of these models is given by

λ t = V N t -(t)h(V N t -(t)), (5)
where h(•) is the failure rate of the first failure time and is assumed to be differentiable. We use the Weibull distribution for the first failure time. Thus, the initial failure rate is a power law process (PLP):

h(t) = β η t η β-1 = αβt β-1 , ( 6
)
where η is the scale parameter and β is the shape parameter. The form with α is another parameterisation that is more convenient in certain computations, with α = η -β . The effect of the maintenance actions on the failure process is characterised by a particular virtual age parameterisation.

Several forms can be assigned to the virtual age. The arithmetic reduction of age model with memory 1 (ARA 1 ) was presented in [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF]. This corresponds to the Type I Kijima model [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] with a deterministic constant maintenance effect. In this case, the maintenance is supposed to reduce the virtual age by an amount that is proportional to the accumulated age since the previous failure, and we obtain

V i (T i ) = V i-1 (T i ) -ρ(V i-1 (T i ) -V i-1 (T i-1 )). ( 7
)
Its virtual age at time t after the ith failure is given by

V i (t) = t -ρT i , (8)
where ρ denotes the efficiency parameter. If ρ = 1, we return to the AGAN case and the maintenance is perfect. If ρ = 0, we are in the ABAO case and maintenance is minimal. The maintenance is efficient between 0 and 1, although not perfect. If ρ is negative, the maintenance is harmful and the system becomes older after maintenance.

For illustration purposes, we use the ARA 1 model in our work, but other virtual age models can be used, such as the ARA ∞ model of [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF] or the quasi-renewal model of [START_REF] Wang | A quasi renewal process and its applications in imperfect maintenance[END_REF], with the parameterisation given in [START_REF] Doyen | A generic framework for generalized virtual age models[END_REF].

Optimal preventive maintenance planning without covariates

Let us first consider the classical case in which the conditional failure intensity is not influenced by any covariate to recall several existing results on which we build our contribution.

Minimal repair and periodic replacement

In the case of minimal repair with periodic preventive replacements, Λ t is no longer random and G(τ ) = C P rev +C Corr Λτ τ . [START_REF] Barlow | Optimum preventive maintenance policies[END_REF] optimised this cost function and showed that τ 0 P = B -1 (C P rev /C Corr ), where B(τ ) := λ τ τ -Λ τ , is the optimal replacement period when the corrective repairs are minimal.

Imperfect repair and dynamic replacement date

Both periodic and dynamic preventive replacement policies can be considered in the presence of imperfect repairs.

Periodic replacement. For non-minimal but imperfect corrective repairs, the failure process between two renewals is no longer an NHPP, and thus, the cumulative intensity is random. Therefore, we let

Φ(t) = E[N t ] = E[Λ t ] and (9)
φ(t) = Φ (t). ( 10
)
We assume that φ(•) is strictly increasing and differentiable for all t ≥ 0. The mean asymptotic cost per time unit for the periodic maintenance policy is now

G(τ ) = C P rev + C Corr Φ(τ ) τ . (11)
Then,

G (τ ) = C Corr φ(τ )τ -(C P rev +C Corr Φ(τ )) τ 2
, where G (τ ) = 0 is equivalent to B(τ

) := φ(τ )τ - Φ(τ ) = C P rev
C Corr (the expression of B(t) defined in Section 2.3.1 corresponds to the particular case in which Λ t is no longer random). By studying the second derivative G (τ ), we can show that the optimal periodic replacement date becomes τ 0

P = B -1 (C P rev /C Corr ), where B(t) := tφ(t) -Φ(t) is monotonic.
Dynamic replacement. In the presence of imperfect corrective repairs, and as illustrated in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], the optimal replacement policy is no longer periodic, but it should use the failure process history. The rationale for an optimal dynamic replacement policy is that if the system fails directly before a preventive replacement scheduled at τ 0 P , efficient corrective maintenance can decrease the conditional system failure intensity, and thus, the preventive replacement should be postponed. To consider the system history, we denote, for all t ≥ 0, the σ-algebra F t = σ{N s , s ≤ t}, and for all 0 ≤ s ≤ t,

Φ(t|F s ) = E[N t |F s ] (12) φ(t|F s ) = lim h→0 1 h (Φ(t + h|F s ) -Φ(t|F s )). ( 13
)
Again, we assume that φ(t|F s ) is strictly increasing and differentiable for all 0 ≤ s ≤ t.

With knowledge of the failure/repair process history until time s, let us consider the preventive maintenance policy until time s + T (T > 0), which consists of n preventive replacements at times s + x 1 , s + x 1 + x 2 , . . . , s + x 1 + . . . + x n with x 1 + . . . + x n ≤ T . [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] demonstrated that the mean conditional cost of such a policy is

C(n; x 1 , . . . , x n |s, T ) = nC P rev + C Corr   Φ(s + x 1 |F s ) -N s + n i=2 Φ(x i ) + Φ   T - n j=1 x j     (14)
and that the optimal policy satisfies

x * 2 = . . . = x * n * = T -x * 1 n * , (15)
φ(s + x * 1 |F s ) = φ T -x * 1 n * , (16)
and

n * = T -x * 1 B -1 (C P rev /C Corr ) + (s, T ), where | (s, T )| < 1. (17) Moreover, for any i ∈ {2, . . . , n * }, lim T →∞ T n * = lim T →∞ x * i = B -1 (C P rev /C Corr ) (18)
and lim

T →∞ φ(s + x * 1 |F s ) = φ(B -1 (C P rev /C Corr )). ( 19
)
This means that after the first preventive replacement (performed at time s + x * 1 ), the following preventive maintenance dates are periodic, with an asymptotic period of τ 0 P = B -1 (C P rev /C Corr ) for an infinite time horizon. This is because each preventive replacement is a renewal, which resets the system history to zero, and after s, no new information will be obtained from the system.

At this point, let us assume that at each time s, F s is known. After each failure/corrective repair event, the previously scheduled time point of preventive maintenance may be revised according to the new information that corrective maintenance has been performed. Using the fact that λ t = φ(t|F t -), it was shown in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] that with this continuous information, the next dynamically updated preventive maintenance date τ 0 D is defined as the first hitting time for the conditional failure intensity:

τ 0 D = inf{t > 0 : λ t ≥ φ(B -1 (C P rev /C Corr ))}, ( 20
)
where the corresponding preventive maintenance intensity threshold is

s λ = φ(B -1 (C P rev /C Corr )). (21)
Subsequently, the dynamic maintenance policy can be implemented in two steps. First, the intensity threshold s λ must be computed. Subsequently, during the system life, the conditional intensity must be continuously evaluated and compared to s λ . As soon as the intensity reaches this threshold, preventive maintenance is triggered. In practice, the conditional intensity is not continuously computed but only determined at each corrective maintenance time point. Indeed, the conditional intensity is a deterministic function of the virtual age, which is only modified by a failure/corrective maintenance event. Therefore, the conditional intensity is dynamically updated and computed with the current virtual age following each corrective maintenance to determine the threshold hitting time, which sets the next preventive maintenance date. This tentative replacement date is then either reached without any other corrective maintenance occurrence (and the system is replaced), or corrective maintenance occurs beforehand and the conditional intensity is computed again to determine the newly planned preventive maintenance date. This dynamic policy is less expensive than the periodic policy corresponding to criterion (11) in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF].

The threshold s λ for a given failure intensity cannot be computed analytically in the general case of imperfect maintenance models. Instead, as proposed in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], it can be estimated by Monte Carlo simulation using the fact that Φ(•) is convex; see Algorithm 1.

Algorithm 1 Computation of preventive maintenance intensity threshold

1: Simulate M failure paths with the same model and the same parameters as the studied system 2: Compute Φ(t) as the greatest convex minorant of the Nelson-Aalen estimator of Φ(t) from the M simulated paths 3: Compute φ(t) as the right derivative of Φ(t) 4: Compute B(t) = t φ(t) -Φ(t); 5: Invert B to compute an estimation of B -1 (C P rev /C Corr ) and then of s λ Fig. 1 depicts the optimisation principle of preventive maintenance dates within this framework. The estimated mean intensity φ(t) is represented by a thin black dotted curve. It reaches the estimated threshold ŝλ , which is denoted by the horizontal grey dotted line, at the time marked by the vertical grey line. This time point corresponds to the theoretical preventive maintenance time τ 0 . However, the observed intensity in red reaches the estimated threshold slightly before τ 0 . Hence, the observed preventive maintenance time instant, which is denoted by the vertical red dotted line, occurs before τ 0 at the time when the observed intensity hits the estimated threshold. After this first preventive maintenance, the observed intensity, denoted in orange, decreases after 15000 time units: a failure occurs and corrective maintenance is performed, thereby lowering the intensity. Therefore, this observed intensity hits the estimated threshold later than the mean intensity, and the second preventive maintenance date, which is denoted by the vertical orange dotted line, is postponed.

Virtual age models with covariates

When designing a dynamic preventive maintenance policy that is adapted to each individual maintained system, it is interesting to consider the differences between the systems owing to their specific failure/repair history as well as the heterogeneity between systems. Such heterogeneity between the systems can be considered with covariates. This information may be constant when it represents, for example, the manufacturer of a system, or dynamic in the case of a covariate corresponding to a degradation process. A modelling framework for virtual age models with timedependent covariates was proposed in [START_REF] Brenière | Virtual age models with time-dependent covariates: A framework for simulation, parametric inference and quality of estimation[END_REF], together with the associated simulation procedure to generate synthetic datasets and an inference method to estimate the model parameters. The hypothesis is that covariates can be represented satisfactorily by a stepwise constant model (in which constant covariates are a particular case). Within this framework, and following a Cox proportional hazard model approach [START_REF] Cox | Regression models and life tables (with discussion)[END_REF], the covariates are included in the conditional failure intensity with the following expression:

λ t = V N t -(t)h(V N t -(t)) exp(γ X t ), ( 22
)
where γ is the coefficient measuring the influence of the covariates and X t , t ≥ 0 is the path of the covariates.

The remainder of this paper presents the contribution of this work, whereby covariate information is introduced into the preventive maintenance decision-making procedure.

Optimal preventive maintenance planning with fixed covariates

In this section, we introduce fixed covariates in the optimisation of the preventive replacement date. Different levels of information on these covariates can be considered: their value may be known, only their distribution may be known, or no information may be available. This information on the covariates may also be available at different stages of the preventive maintenance decision-making procedure (threshold computation or intensity evaluation). Thus, we study how the information on the covariates can be integrated into the decision and lead to different replacement times depending on the available information, as well as the impact of the information level on the cost of the preventive replacement policy.

Minimal repair

We first consider the case of minimal repair, and assume that for each system, the covariates take a new value after each preventive maintenance (renewal). The failure intensity between the two preventive renewals is expressed as

λ t = exp(γ X)h(t) (23)
and the cumulative intensity is

Λ t = exp(γ X)H(t), ( 24
)
where X is the covariate value after the last renewal, γ is the coefficient of influence of the covariates, h(•) is the failure rate of the first failure time without covariates, H(•) is its primitive, which is null in 0, and t is the time elapsed since the last renewal. Two cases must be distinguished depending on the available level of information on the covariates at the beginning of the system life; that is, at the time of computation of the maintenance threshold. The value χ that is taken by the covariates after the different successive renewals; that is, at the beginning of each cycle, can either be known or not.

Covariates as a priori unknown random variables

In the case in which the value χ of the covariates is unknown at the time of the threshold computation, only the knowledge of its probability law can be used in the failure intensity expression. Thus, we obtain

Φ(t) = E[Λ t ] = E[exp(γ X)]H(t) (25)
and

φ(t) = Φ (t) = E[exp(γ X)]h(t), (26)
because H(•) is a deterministic time function. We can then write

B(t) = tφ(t) -Φ(t) = E[exp(γ X)](th(t) -H(t)). ( 27
)
By writing B * (t) := th(t) -H(t), we obtain

B -1 (t) = B -1 * t E[exp(γ X)] , (28)
and thus, the maintenance threshold expression for the failure intensity becomes

s λ = φ(B -1 (C P rev /C Corr )) = E[exp(γ X)]h B -1 * C P rev /C Corr E[exp(γ X)] . ( 29
)
The observed failure intensity can now be compared to this threshold s λ to determine the preventive maintenance date, and two cases must be distinguished depending on whether the covariate values are known at the decision-making time; that is, at the time of this comparison.

Covariate value unknown at decision-making stage. If the new value χ of the covariate is still unknown when the system is returned to work after preventive replacement, the expected value of the failure intensity (innovation theorem; see [START_REF] Andersen | Statistical models based on counting processes[END_REF]) must be used to determine the preventive maintenance date τ f ix

P such that E[exp(γ X)]h(τ f ix P ) = E[exp(γ X)]h B -1 * C P rev /C Corr E[exp(γ X)] , (30)
and hence,

τ f ix P = B -1 * C P rev /C Corr E[exp(γ X)] , (31)
because we assume that h(•) increases. Note that τ f ix P is periodic, as nothing modifies its value during the life of the system. This result can be compared to the case in which the knowledge of the covariate distribution is not considered; that is, when it is considered that γ = 0, which leads to

τ f ix P = B -1 * (C P rev /C Corr ). ( 32
)
τ f ix P is a periodic maintenance time point, similar to τ f ix P ; however, compared to expression (32) of τ f ix P , expression (31) of τ f ix P clearly shows the effect of the covariates on the preventive maintenance date.

Covariate value known at decision-making stage. Conversely, if the value χ of the covariates after each preventive maintenance (renewal) is known, it can be integrated into the failure intensity expression, and the preventive maintenance date τ f ix D is such that

exp(γ χ)h(τ f ix D ) = E[exp(γ X)]h B -1 * C P rev /C Corr E[exp(γ X)] , (33)
and hence,

τ f ix D = h -1 E[exp(γ X)] exp(γ χ) h B -1 * C P rev /C Corr E[exp(γ X)] . (34)
In the case with covariates, the structure of this maintenance policy corresponds to that proposed by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] and is summarised in subsection 2.3.2. Indeed, the obtained policy is dynamic and no longer periodic because the covariate value is known only after preventive maintenance; that is, after system renewal. In this case, again, comparing expression (34) of τ f ix D with expression (20) of τ 0 D clearly demonstrates the effect of the available knowledge on the covariates on the preventive maintenance date determination.

Covariates as a priori known random variables

In a full-information situation, the value χ is known at both the threshold computation time and decision-making time. By using this value in the failure intensity expression, we obtain

Φ(t) = exp(γ χ)H(t) (35) and φ
(t) = exp(γ χ)h(t). (36)
Hence,

B(t) = exp(γ χ)B * (t), (37)
and thus,

sλ = φ( B-1 (C P rev /C Corr )) = exp(γ χ)h B -1 * C P rev /C Corr exp(γ χ) . ( 38
)
By comparing the observed intensity (thus, with χ) to this threshold, we obtain the expression of the preventive maintenance date for this fully informed policy:

τ f ix D = B -1 * C P rev /C Corr exp(γ χ) . ( 39
)
As the covariate value χ changes at each renewal, the replacement policy τ f ix D is not periodic.

Comparison of four replacement policies

A numerical comparison of the maintenance costs that are incurred by the four replacement policies is presented in this section. For comparison purposes, Table 1 displays the covariate knowledge levels for each policy.

Experimental protocol. To compare the mean cost per time unit incurred by the four maintenance policies, N system histories were simulated until their first preventive maintenance dates, which was equivalent to simulating a single system undergoing N preventive replacements, as they were supposed to be perfect. Following (1), the total cost of system j = 1, . . . , N until its first observed preventive maintenance date τ (j) was equal to where N (j) τ (j) denotes the number of observed failures in system j until time τ (j) . Following (2), the mean asymptotic cost per time unit was estimated from the simulation results of the N system histories, as follows:

C (j) (τ (j) ) = C P rev + C Corr N (j) τ (j) , (40)
τ f ix P τ f ix P τ f ix D τ f ix D E[exp(γ X)]
Ĝ(τ ) = N j=1 C (j) (τ (j) ) N j=1 τ (j)
.

(41)

This procedure was repeated for the four policies with N = 100 and a Weibull initial intensity of η = 15000 for the scale parameter and β = 2 for the shape parameter. With such an initial intensity, B -1 * (t) could be explicitly computed and it was not necessary to estimate Φ(t) based on simulations, as in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF]. The preventive renewal and minimal repair costs were C P rev = 1 and C Corr = 3 (in arbitrary units), respectively.

A single covariate was considered, for which the value on a system between two successive renewals followed a uniform distribution on {-5, -4, . . . , 4, 5}, which means that for each system history, a new covariate value was drawn from this uniform distribution in the simulation. Of course, this covariate value was used differently according to the different policies and the associated assumed level of covariate knowledge. Furthermore, to study the influence of coefficient γ on the maintenance cost, different values of γ ∈ {-2, -1, 0, 1, 2} were considered.

Numerical results and discussion. of γ arose from the fact that the possible covariate values were symmetric and had a uniform distribution on {-5, -4, . . . , 4, 5}. Thus, taking γ = 1 or γ = -1 amounted to taking the opposite of the covariate, and because the same weight existed on both the negative and positive covariate values, a symmetry appeared in the costs. This observation of symmetry is relevant to all policies. Moreover, when γ > 0, the policy that yielded τ f ix P (red curve) performed a replacement too late (too many failures would occur beforehand) if χ > 0 and too early (the preventive maintenance would be too hasty) if χ < 0. This was reversed if γ < 0. In both cases, the cost was higher than when the covariate had a null effect (i.e. γ = 0). This behaviour was observed for all maintenance policies except for τ f ix D , which is discussed in Appendix A. Moreover, the mean costs per time unit in all policies were equal for γ = 0, which can be explained by examining the expression of the preventive maintenance date in each policy and replacing γ with 0, and it can be observed that τ

f ix P = τ f ix P = τ f ix D = τ f ix D = B -1 * (C P rev /C Corr ).
Apart from γ = 0, the curves piled up into a clear classification. The periodic policy, which considered E[exp(γ X)], thereby yielding τ f ix P (orange curve), was better than that of [START_REF] Barlow | Optimum preventive maintenance policies[END_REF], resulting in τ f ix P (red curve), which considered no covariate information, especially when the absolute value of γ was large. Thus, the integration of the knowledge of the covariate distribution into the maintenance decision had a positive impact on the costs. Furthermore, the dynamic policy corresponding to the approach of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] adapted to covariates, resulting in τ f ix D (green curve), was better than the previous one, where the observed value of χ was not considered at all. In this case, the optimality theorem of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] proves that the policy τ f ix D was superior in terms of cost compared to τ f ix P . Therefore, observing χ is useful for planning preventive maintenance. Finally, the blue cost curve, which corresponded to τ f ix D , and therefore, to the case in which χ was known since the threshold computation stage, was the lowest, demonstrating that this policy was the best and incurred the lowest cost. These observations demonstrate that integrating at best the monitoring information that is available online (in this static case, the covariate value after a preventive renewal) in the maintenance decision-making enables the optimisation of the preventive maintenance replacement time points. [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] proposed a dynamic policy adapted to imperfect corrective repairs but without considering covariates. It can be observed in this case that when the covariates intervened, this policy was no longer the best and our proposed policy outperformed it in the presence of covariates.

Imperfect repair

Now, consider the case in which corrective maintenance actions are implemented by imperfect repairs. In this case, the observed system failure times also provide valuable information on the failure intensity evolution, and hence, for maintenance decision-making. Let us assume that these dates are known progressively for decision-making during the system life. In this setting, closer to that adopted in Gilardoni et al. ( 2016) but with a system influenced by fixed covariates, we compare the cases in which the covariates are completely unknown (which yields the τ f ix P policy) and the case in which they are known at the threshold computation stage (which yields the τ f ix D policy).

Maintenance decision-making. In this new setting including imperfect repairs, we denote

B * (t) = tE[V N t -(t) h(V N t -(t))] -E[ t 0 V N s -(s) h(V N s -(s))ds]. ( 42
)
When the covariates are observed, B(t) = exp(γ χ)B * (t); when they are not observed, B(t) = B * (t). The computation of B -1 * (t) is no longer explicit and Φ *

(t) = E[ t 0 V N s -(s)h(V N s -(s)
)ds], which is the expected value of the cumulative intensity without covariates, is estimated using the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF]. For the policy with unknown covariates, the preventive replacement date is determined as the time point at which the conditional intensity without covariates reaches the threshold φ * (B -1 * (C P rev /C Corr )). For the policy in which the covariates are known, exp(γ χ) is a factor of both the failure intensity and its threshold, which is simply computed as the time point at which

V N t -(t)h(V N t -(t)) reaches the threshold φ * B -1 * C P rev /C Corr exp(γ χ)
.

Numerical study. For the numerical study, we considered the ARA 1 model with ρ = 0.5 as the repair efficiency parameter2 and a Weibull initial intensity with η = 15000 and β = 2. The maintenance costs were C P rev = 1 and C Corr = 3 (in arbitrary units), the coefficient γ ∈ {-2, -1, 0, 1, 2}, and the covariate took uniformly distributed values in {-5, -4, . . . , 4, 5}. A total of N = 100 system histories were simulated. Fig. 3 presents the mean costs per time unit obtained in the two cases of covariate knowledge with respect to γ. Symmetry can be observed in the cost curves and the replacement policies were equivalent for γ = 0. Furthermore, the information contributed by the covariate was crucial, as it enabled a drastic reduction in the maintenance cost if it was properly integrated into the maintenance decision process; the mean cost per time unit was approximately 17 times higher when γ = 1 and more than 66500 times higher when γ = 2.

Conclusions

In this section, we demonstrated that the knowledge of fixed covariates that have an influence on a system enables the development of preventive maintenance policies that incur lower costs when the covariates are known and properly used in the decision-making procedure. Therefore, the information contained in the covariates is of high value for the optimisation of preventive maintenance replacements that are performed on a system, and must be used at best, as progressively available, in the maintenance decision process.

Optimal preventive maintenance planning with dynamic covariates, minimal repair

In this section, we again consider the optimisation of preventive maintenance dates, but for time-dependent dynamic covariates corresponding to degradation or wear indicators. The dynamic information contained in the monitored covariates is subsequently considered as it becomes available to plan preventive replacements, following a condition-based maintenance approach. This section also assumes minimal repair for corrective maintenance actions, whereas Section 5 deals with imperfect repair.

Preventive maintenance decision with dynamic covariates

Assumptions and decision principle

We assume that the entire covariate path is not observed, but that the covariates are monitored through periodic inspections. Moreover, we assume that the covariates are not observed upon a failure. At each inspection, the observed covariate value is used to compute and update the failure intensity in the same manner as in the failure process history in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF]. If the observed intensity reaches the threshold s λ before the next inspection, preventive maintenance is triggered; otherwise, no action is performed on the system and the next inspection is scheduled. This maintenance decision process extends that proposed in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] (where after each failure, the virtual age and the conditional intensity are updated to determine preventive maintenance) to integrate the information from the covariate monitoring. To this end, we also assume that the covariates are positive and monotonically increasing, that they are restored to zero by a preventive replacement (as is indeed the case if the covariate corresponds to, e.g. a degradation or wear indicator) and subsequently maintain the same distribution, and that the paths between two preventive replacements are independent. In [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], an analytical proof that the proposed dynamic maintenance policy is better than the periodic policy was presented. This proof relies on a property of the conditional intensity that cannot be proven in the case of the dynamic covariates considered in the present work. Therefore, the performance of our proposed policy was assessed through numerical simulations, and it was also compared to the performance of a periodic policy corresponding to the policy obtained by using our maintenance optimisation procedure without any covariate inspection; that is, with only the knowledge of the covariate distribution. The first replacement time instant was determined from the beginning of the system life, and because the preventive maintenance was a renewal, the replacement dates occurred at fixed intervals.

The principle of the preventive maintenance optimisation based on knowledge of the covariate value is presented in Fig. 4 for a system with minimal repairs. The mean intensity, denoted by the thin black dotted curve, reached the threshold, denoted by the grey dotted line, at the theoretical preventive maintenance date τ 0 , denoted by the vertical grey line. However, the observed intensity, indicated in red, only reached the threshold later. The effective preventive maintenance date, denoted by the vertical red line, was thus postponed owing to the knowledge of the true covariate path. Conversely, after the first preventive renewal, the observed intensity, indicated in orange, increased faster than the mean intensity. The preventive maintenance date, denoted by the vertical orange line, therefore occurred sooner than initially scheduled.

Integration of covariates into decision

To illustrate the integration of the covariates into the decision process, let us assume a single covariate following a stochastic Gamma process with shape and rate parameters a and b, respectively, sampled and held at a period dt cov , which means that the covariate values at kdt cov are the values of a Gamma process on these points, and maintained constant in the interval [kdt cov ; (k + 1)dt cov [. In practice, considering such a piecewise constant covariate is not restrictive (because dt cov can be selected to be arbitrarily small), and it enables the use of the failure simulation algorithm proposed in [START_REF] Brenière | Virtual age models with time-dependent covariates: A framework for simulation, parametric inference and quality of estimation[END_REF], which only works with stepwise constant covariates.

Let X t be the covariate value at time t. The system failure intensity between the two preventive renewals is expressed as

λ t = h(t) exp(γX t ), ( 43
)
where γ is the coefficient measuring the covariate influence and t is the time elapsed since the last renewal.

Recall that l (t) := t/dt cov is a positive integer such that

(l (t) -1)dt cov < t ≤ l (t) dt cov . ( 44
)
At time t = 0, the covariate value is assumed to be X 0 = 0 and it remains the same until t = dt cov . After inspection at I 1 > 0, the covariate value X I 1 is known for the entire time interval [(l (I 1 ) -1)dt cov , l (I 1 ) dt cov ]. However, in the interval ]0, (l (I 1 ) -1)dt cov [, the covariate path remains unknown. At this stage, there are two possibilities to predict this path: either assuming that the covariate is constant between two inspections and takes the last known value on the interval, or using the knowledge on the covariate distribution. In the second case, we assume that the observed covariate has the following value:

X t = X I i + a b (l (t) -l (I i ) )dt cov (45) for all t ∈ [(l (I i ) -1)dt cov , (l (I i+1 ) -1)dt cov [
, where I i and I i+1 are two successive inspection time points, with dt cov I i+1 -I i for all i ≥ 1. This corresponds to taking the mean covariate path from I i , knowing that the covariate had the value X I i at time I i and that it moved step by step, as illustrated in Fig. 5. We adopted the second approach in the following study.

Computation of threshold on conditional failure intensity

In the case of minimal repairs, as considered in this section, and in the presence of a covariate, two methods are available for computing the intensity threshold s λ (see Appendix B): one proposed by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] ("Method 1") and our proposition ("Method 2"), which consists of estimating φ(t) and Φ(t) directly from the simulated covariate paths and then computing B(t) = t φ(t) -Φ(t), from which the threshold s λ on the conditional failure intensity can be estimated as

(t) l (Ii) )dtcov ⌧l (I i ) 1 Ii ⌧l (I i ) ⌧l (t) 1 t ⌧l (t) ⌧l (I i+1 ) 1 Ii+1 ⌧l (I i+1 ) = = = = = = (l (Ii) 1)dtcov l (Ii) dtcov (l (t) 1)dtcov l (t) dtcov (l (Ii+1) 1)dtcov l (Ii+1) dtcov dtcov XI i XI i+1 1
Figure 5: Inspection dates and assumed covariate value ŝλ = φ( B-1 (C P rev /C Corr )). Both computation procedures of this threshold allow the integration of information on the dynamic covariate behaviour, which enables the replacement date to be determined dynamically. The results using both methods are reported and compared in the numerical experiments section.

Dynamic determination of replacement date

The replacement date is determined differently depending on whether or not the dynamic covariate value is known by inspection.

• If there is no covariate inspection, it is supposed to follow the mean path provided by its distribution, and the preventive maintenance date is periodic and equal to

τ dyn P = B -1 (C P rev /C Corr ). (46)
• In contrast, if the covariate is known by inspection with time step dt ins , the knowledge of the real covariate value (and no longer of its mean path only) at an inspection point enables the re-evaluation and adaptation of τ dyn D upwards if the observed value is lower than the mean value and downwards otherwise (in both cases with a positive γ).

Knowing that the covariate level at the last inspection time point t I is χ and that the next inspection will take place at time t next I , the next preventive maintenance date can be (re-)planned after each inspection following Algorithm 2.

Algorithm 2 Preventive maintenance date optimisation on system with dynamic covariates, minimal repair 1: The current covariate step is l ← t I /dt cov and it will be l next ← t next I /dt cov at the next inspection

2: if h(t I ) exp(γχ) ≥ s λ then 3:
The next preventive maintenance is planned immediately

4: else if h(t next I ) exp(γ(χ + (l next -l) * dt cov * a/b)) ≥ s λ then 5: l P M ← l+1 6: while h(dt cov * l P M ) exp(γ(χ + (l P M -l) * dt cov * a/b)) < s λ 7: l P M ← l P M +1 8:
end while

9:

The next preventive maintenance is planned at time h -1 (s λ / exp(γ(χ+(l P M -l) * dt cov * a/b))) 10: else 11:

No preventive maintenance is planned before the next inspection 12: end if

Maintenance optimisation: numerical experiments

The proposed dynamic maintenance policy was evaluated through numerical experiments, and two means of optimising its performance were investigated: i) the optimisation of the covariate inspection period and ii) the optimal scheduling of a single inspection between two replacements.

Optimal covariate inspection period

Experimental protocol. In this first numerical experiment for studying the effect of the inspection period on the policy performance, the inspection period was selected with respect to τ 0 = B -1 (C P rev /C Corr ), the mean time for the intensity to hit threshold s λ , with the aim of observing a number of inspections close to the a priori selected value. This procedure enabled the policy performance assessment on a set of simulated histories that were consistent in terms of the number of performed inspections. If one inspection was sought during a maintenance cycle (i.e. between two replacements), the inspection period was taken as dt ins = τ 0 /2, and in most of the systemsimulated histories, one inspection occurred around the middle of the system life. Other inspections could occur if τ dyn D exceeded τ 0 . If n inspections were sought on average, the inspection period was taken as dt ins = τ 0 /(n + 1). In the numerical study, cases with n ∈ {1, 2, 4, 8, 16}, in addition to the case without inspections (i.e. n = 0), were considered.

Without inspection cost. We first consider the case in which inspections are performed at no cost; it is advantageous to inspect as often as possible, so as to follow the covariate evolution closely and update the conditional failure intensity to schedule the best preventive maintenance date, thereby avoiding too many failures or too early preventive maintenance, which would increase the total cost in both cases. This situation was simulated with η = 15000 and β ∈ {1.5, 2, 2.5, 3} for the PLP parameters, with different maintenance cost settings (C P rev = 1 and C Corr ∈ {3, 5, 15} or C P rev ∈ {3, 5, 15} and C Corr = 1, all in arbitrary units), γ = 1 for the covariate influence, a = 10 -3 , b = 8, and dt cov = 10 for the covariate simulation, M = 10000 for the number of simulations for the threshold computation, n ∈ {0, 1, 2, 4, 8, 16} for the inspection number, and N = 10 6 for the number of simulated system histories for estimating the mean total cost. It was verified that with these parameters, dt insp = τ 0 /(n + 1) lay between 234 and 2725 time units, which was still much larger than the covariate discretisation step dt cov = 10. Fig. 6a presents a characteristic example of the mean cost per time unit that was obtained in one of the cost configurations with respect to the number of inspections and with the two methods of threshold computation. The red curve corresponds to the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], whereas the blue curve corresponds to our method. As expected, a lower inspection period (i.e. more inspections on the system life) resulted in a lower maintenance cost, which enabled the quantification of the covariate information impact on the system preventive maintenance date optimisation. However, not all of the studied configurations had exactly the same type of cost curves; in some cases, instabilities occurred (i.e. the cost curve was not smoothly decreasing), even with N = 10 6 simulated system histories, which is visible from the blue curve in Fig. 6a. It was found that such instability was heavily dependent on the parameters, and the curve in Fig. 6a represents a particular case in which the cost was smoothly decreasing. For example, a lower ratio C P rev /C Corr led to more unstable results. In fact, if C Corr (the corrective maintenance cost) was high, preventive maintenance had to be performed as early as possible to avoid failures. However, the inspection dates were set with respect to the theoretical preventive maintenance date τ 0 , and if n inspections were sought before τ 0 , they were planned periodically with period τ 0 /(n + 1). Thus, if τ 0 was low, the inspections were close to one another, even when n was not high, and therefore, there was no striking difference between the different inspection periods and the curve cost would possibly not be smoothly decreasing.

Regarding the threshold computation method, in most cases, our proposed method led to lower costs than the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] (see Appendix B), although it remained unclear in certain configurations owing to the instabilities.

With inspection cost. We consider the case in which the inspections have a cost C I . It would be interesting to optimise the inspection period to determine the optimal trade-off between too frequent inspections (which would be too expensive with respect to the gain resulting from the additional information) and too sparse inspections (which would not allow proper adaptation of the maintenance decision). The simulation study presented above was repeated, this time including the inspection cost. The total cost between the two preventive replacements was

C(τ ) = C P rev + C Corr N τ + C I N (I) τ , (47)
where

N (I) τ
denotes the number of inspections performed during the cycle. Fig. 6b presents the cost curves that were obtained with the same configuration as in Fig. 6a, but with an inspection cost set to C I = 0.003 (in arbitrary units). These cost curves clearly indicate the minimum for n = 2, which was the optimal number of inspections before τ 0 in this configuration. Beyond this number, new inspections would be too costly with respect to the value of the information they contributed, and below this number, information on the system would not be sufficient.

Figs. 6a and 6b illustrate the fact that the proposed preventive replacement decision rule enabled the efficient use of the available covariate information to schedule the preventive replacement date, leading to lower maintenance costs than a covariate-blind decision rule (i.e. with n = 0 inspection).

Covariate with random effect. If the cost curve shown in Fig. 6b exhibited a clear minimum, the improvement was only 2.72% between the cost in the case without inspection (i.e. a periodic replacement policy) and the minimum cost reached with two inspections. This low improvement in the cost was assumed to be owing to the too low variability of the covariate from one path to another; in such a situation, monitoring the covariate returned information with limited added value. To increase the variability in the covariate, a new simulation study was established, in which the covariate was modified to add more path-to-path variability. Each path started from a random centred value, with a random effect on the shape parameter a. The computation of the threshold used both our method and that adapted from [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], with a mean value of a and assuming no offset at the origin. Note that this mean value of a was the only data available on this modification of the covariate; during the optimisation of the preventive replacements, covariate path prediction was still performed with the observed value of the covariate and the mean value of a. However, each simulated path of the covariate used a value of a and an offset at the origin that was randomly selected before the beginning of the system life.

The simulation was performed with β = 2, C P rev = 1, C Corr = 3 (in arbitrary units). The offset was uniformly drawn between {-1, -0.5, 0, 0.5, 1} and a between {0.001, 0.002, . . . , 0.005}. more marked, flattening the remainder of the curve. A magnified view of this curve is shown in Fig. 7b. The cost improvement on the blue curve (obtained using our proposed threshold computation method) between the case without inspection (i.e. periodic replacement policy) and that with two inspections was approximately 21.0%. This higher cost improvement demonstrates the fact that in a situation of higher variability in the covariate behaviour, covariate inspection contributes information with greater added value for maintenance optimisation purposes.

Optimal covariate inspection date

Consider the issue of optimally scheduling a single inspection on the entire system life; that is, before its replacement. Intuitively, an inspection at the beginning of the system life does not have the same impact as when the system has already worked for some time. Initially, the covariate may not have evolved along its path to offer useful information. Conversely, if the inspection occurs too late, the covariate may have exceeded a critical level and preventive maintenance would be scheduled too late, after several costly failures of the system. The issue of determining the best inspection scheduling time was investigated through Monte Carlo simulations.

The simulation study was run with the same parameters as in Section 4.2.1, devoted to the optimisation of the inspection period. It was again assumed that the corrective maintenance action was more expensive than the preventive maintenance action (C Corr > C P rev ). The possible inspection dates were defined in proportion to the theoretical preventive maintenance date τ 0 ; the inspection could occur in 0.1τ 0 , 0.2τ 0 , . . . , τ 0 . The algorithm that was used to deal with this situation was the same as Algorithm 2, with lines 4 (but not 5 to 9) and 10 to 12 removed because there would not be a next inspection.

As a characteristic configuration, Fig. 8 presents the mean cost per time unit as a function of the date of inspection (more precisely, of the proportion of τ 0 ), and with the two threshold computation methods. The cost curve exhibits a minimum that was observed for an inspection at 0.8τ 0 using our threshold computation method.

It should be noted again that the minimum is not clearly marked for all parameter configurations, and unsmooth cost curves could be obtained when C P rev /C Corr was low. To explain this phenomenon, recall that the date of a single inspection was set with respect to the theoretical preventive maintenance date τ 0 : the inspection was performed at x % of τ 0 , with x ∈ {10, 20, . . . , 100} (the different values of x constitute the cost curve abscissa). Hence, if C Corr was high, τ 0 was low, and the covariate variations were not sufficiently significant between the different values of x to constitute valuable information, which led to a less marked minimum.

Optimal preventive maintenance planning with dynamic covariates, imperfect repair

This section considers the optimisation of preventive maintenance dates for time-dependent dynamic covariates and imperfect repairs, which is the most complete situation investigated in this study. Imperfect repairs are assumed to have no effect on the covariate's level and evolution: in contrast to replacements, they modify the failure intensity without affecting the covariate itself.

Preventive maintenance decision procedure

In the case of imperfect repair, our proposed method for threshold computation is no longer relevant and the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], including covariate simulations in the compu-tations, must be used. As this method is more computationally demanding, the simulations were performed with a larger covariate simulation time step, dt cov = 100.

Algorithm 3 outlines the modified procedure to optimise the preventive maintenance dates in the case of imperfect repair and dynamic covariates. As an imperfect repair modifies the conditional failure intensity, and contrary to Algorithm 2, in this case, the preventive maintenance date is modified at each corrective repair, even though the covariate is not observed. The remainder of the algorithm is very similar. It is presented for an ARA 1 virtual age model, but it can easily be adapted to other virtual age models. Knowing that

• the covariate level at the last inspection time point t I is χ,

• the next inspection will take place at time t next I , and

• the last corrective maintenance took place at time T CM , the next preventive maintenance date can be (re-)planned after each inspection or corrective maintenance following Algorithm 3.

Algorithm 3 Preventive maintenance date optimisation on system with dynamic covariates, imperfect repair 1: The current covariate step is l ← t I /dt cov and it will be l next ← t next I /dt cov at the next inspection

2: The current time is t ← max(t I , T CM ) 3: if h(t -ρT CM ) exp(γχ) ≥ s λ then 4:
The next preventive maintenance is planned immediately No preventive maintenance is planned before the next inspection 16: end if

Optimal covariate inspection period: numerical experiment

In this section, as in Section 4.2.1, the inspection period optimisation is investigated. The inspection period dt ins was set with respect to τ 0 = B -1 (C P rev /C Corr ); i.e., dt ins = τ 0 /(n + 1), with n ∈ {1, 2, 4, 8, 16}, in addition to the case without any inspection.

Without inspection cost

Consider the case with no inspection cost. The numerical experiments were performed with same parameters as in Section 4.2.1, but with an ARA 1 imperfect maintenance model with an efficiency parameter ρ ∈ {0.1, 0.3, 0.5, 0.7, 0.9}. As discussed previously, only the case in which C Corr > C P rev (repair is more expensive than preventive replacement) was considered. Fig. 9a presents the mean cost per time unit as a function of the number of inspections within a replacement cycle for a characteristic configuration example. It can be observed that the maintenance cost rate decreased with the number of inspections. Some parameter configurations were more stable than others, but the global curve shape remained the same in all studied configurations. This work is at the confluence of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] and Section 4.2.1, and the results show how both the covariate information and failure/repair history knowledge can be used with profit in an imperfect maintenance context with dynamic covariates to schedule preventive replacement more effectively.

With inspection cost

Consider the case in which inspections are performed at a cost C I , as in Section 4.2.1. Fig. 9b presents the costs that were obtained in the same configuration as in Fig. 9a, but with an inspection cost set to C I = 0.005 (in arbitrary units). The cost curve clearly shows a minimum that corresponds to the optimal inspection period for the system considered. With fewer inspections, important information on the covariate evolution would be missed. More inspections would result in an excessive total inspection cost with respect to the gain in terms of maintenance costs. In particular, Fig. 9b compares the proposed dynamic replacement policy integrating the covariate information with a covariate-blind one proposed by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] (obtained for the n = 0 inspection case).

Conclusions

A procedure to optimise preventive replacements in a system with fixed or dynamic covariates and with minimal or imperfect repairs has been proposed. This procedure is implemented in two steps. First, the knowledge of the statistical properties of the covariate (or its value if it is fixed) is used to compute an intensity threshold, beyond which preventive maintenance is to be triggered. Subsequently, during the system life, the conditional intensity, which is computed with the true values of the observed covariates (known at some inspection points), is compared to the threshold, and preventive maintenance is performed as soon as it is exceeded.

The numerical experiments demonstrated that the knowledge of the covariate enabled the cost performance of the preventive replacement policy to be improved compared to the periodic "covariate-blind" policies (for minimal repairs, [START_REF] Barlow | Optimum preventive maintenance policies[END_REF]) or dynamic ones (for imperfect repairs, [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF]). The influence of the level of covariate information on the cost-optimal solution was also highlighted. Therefore, the use of dynamic information on failures and covariates is of great interest to optimise preventive replacement scheduling and to reduce maintenance costs. In particular, for costly inspections, it was found that it is possible to determine an optimal covariate inspection period or the optimal time to schedule an inspection when a single inspection can be performed during the system life.

This final point opens the research field. Indeed, if it is possible to implement a single inspection, a second, third, and so on could be performed until it no longer improves the total cost; that is, when too many inspections are performed. This could allow us to depart from periodic inspection policies and determine the optimal number of inspections to be performed along the system life, as well as their dates.

Furthermore, in this study, we assumed that the virtual age model parameters, those of the initial intensity, and the covariate coefficient were known. However, the estimation of these parameters would introduce an additional bias into diverse optimisations, which would be interesting to investigate in a deeper analysis. Moreover, the parameters of the covariate stochastic process could be estimated dynamically based on the information provided by inspections to develop a dynamic adaptive maintenance optimisation procedure [START_REF] Liu | A dynamic prescriptive maintenance model considering system aging and degradation[END_REF]). higher γ results in a worse effect on the system. We observe that for a centred range of covariates, the curve lies in between but is closer to that of the positive range when γ is negative and closer to that of the negative range when γ is positive, which results in the overall appearance of the plain curve. Furthermore, we can examine the mathematical expression for the mean cost per time unit. With a Weibull distribution for the first corrective maintenance time point, as in this study, the expected value of the number of failures that occurred before the preventive maintenance time point is constant and depends on neither γ nor the covariate value. Indeed, because the corrective repairs are minimal, where Λ(t) = t 0 λ t denotes the cumulative intensity of the failure process. Then, the numerator of the cost in (3) is constant. The denominator is simply the expected value of the preventive maintenance date, the expression of which is presented in (39). range case (plain curve) lies between the two other cases. For a positive covariate, the curve decreases: a higher γ resulted in a smaller preventive maintenance time point because the system was wearing faster. This is the opposite of the negative covariate. These elements may explain the non-convex blue curves plotted in Fig. 2 Thus, two methods are available for estimating the intensity threshold s λ in the case of minimal repair: the general method proposed by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] (see Algorithm 1) and the proposed method. We can compare them from different perspectives: the computation time, value of the threshold, or costs achieved with such thresholds. We now provide insights into the first two perspectives, whereas the costs are examined in Section 4.2.1. For both methods, we computed the same threshold 100 times with η = 15000, β = 2, C P rev = 1, C Corr = 4 (in arbitrary units), dt cov = 100, a = 0.001, and b = 8. Fig. B.12 presents boxplots of the threshold obtained using the two methods, where 1 (red box) represents the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] and 2 (blue box) indicates our new method. Moreover, Fig. B.12 presents boxplots of the computation times using the two methods. First, we observe that the dispersion was much higher in the method of .12: Intensity thresholds and computation times of threshold with two different methods: proposed method in blue and Gilardoni's method in red [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] than in our method. This is because of the simulation of the failure times.

In our experiment, we simulated M = 10000 failure paths. We could have used more time, but the computation would have been longer. Furthermore, we observe that the medians were not equal; they were 8.82 × 10 -5 in the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] and 8.49 × 10 -5 in our method. However, our method remained in the lower whisker of the red boxplot, which is one time of the interquartile range. Regarding the computation time, the results are clear: the method of [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] was by far the longest, with a median of 9.56 s, whereas the median in our method was 4.78 × 10 -1 s. Therefore, our method was 20 times faster.

Figure 1 :

 1 Figure 1: Schematic optimisation principle of preventive maintenance dates with imperfect repairs

Figure 2 :

 2 Figure 2: Mean asymptotic cost per time unit with respect to used policy and γ, minimal repair

Figure 3 :

 3 Figure 3: Mean asymptotic cost per time unit with respect to used policy (known vs. unknown covariate) and γ, imperfect repair

Figure 4 :

 4 Figure 4: Schematic of optimization principle for preventive maintenance dates with covariates

Figure 6 :

 6 Figure 6: Example of cost per time unit, β = 1.5, CCorr = 3, CP rev = 1, minimal repair (vs. number of inspections)

  Figs. 7a and 7b show the mean cost per time unit with the inspection cost C I = 0.003 as a function of n. The cost improvement between the cases with and without inspection was much

Figure 7 :

 7 Figure 7: Cost per time unit with inspection cost, β = 2, CCorr = 3, CP rev = 1, minimal repair, covariate with random effect (vs. number of inspections)

Figure 8 :

 8 Figure 8: Example of cost per time unit, β = 2, CCorr = 3, CP rev = 1, minimal repair (vs. date of inspection)

5:

  else if h(t next I -ρT CM ) exp(γ(χ + (l nextl) * dt cov * a/b)) ≥ s λ then 6: l P M ← l+1 7: while h(dt cov * l P M -ρT CM ) exp(γ(χ + (l P Ml) * dt cov * a/b)) < s λ The next preventive maintenance is planned at date h -1 (s λ / exp(γ(χ + (l P Ml) * dt cov * a/b))) + ρT CM 11:if a failure occurs before the planned preventive maintenance then

Figure 9 :

 9 Figure 9: Example of cost per time unit, ρ = 0.5, β = 2, CCorr = 3, CP rev = 1, imperfect repair (vs. number of inspections)

Figure A. 10 :

 10 Figure A.10: Mean asymptotic cost per time unit with respect to covariate range and to γ, minimal repair

B

  * (t) = th(t) -H(t) = αβt β-1αt β = α(β -1)t β ;

  Fig. A.11 plots the mean preventive maintenance time instant for the three covariate range cases. We again observe that the centred

Figure A. 11 :

 11 Figure A.11: Mean preventive maintenance time point with respect to covariate range and to γ, minimal repair

  l 0 -1)dt cov ) β g l 0 + l 0 -1 l=1 g l ((ldt dov ) β -((l 0 -1)dt cov ) β )

  Figure B.12: Intensity thresholds and computation times of threshold with two different methods: proposed method in blue and Gilardoni's method in red

Table 1 :

 1 Level of knowledge of covariates for computation of different preventive replacement time instants

	known	no yes yes yes
	χ known for threshold computation no	no	no yes
	χ known for decision-making	no	no yes yes

We draw the attention of the reader to the fact that the θ parameter of[START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] corresponds to 1 -ρ

Acknowledgments

This work was partially supported by LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01), funded by the French programme Investissement d'avenir. Most of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir, supervised by the Agence Nationale pour la Recherche.

Appendix A. Elements on shape of curve cost as function of γ for fixed covariate -Fig. 2 In this Appendix, details are presented regarding the case with known fixed covariates and minimal repair, as discussed in Section 3.1.2 and studied (with other policies) in Section 3.1.3.

To understand the reason for the blue curve in Fig. 2 (the cost rate curve as a function of γ in the case of a completely known covariate, i.e. when the replacement date τD ) is not convex, we replicated the study but with different ranges of covariate values: we selected either a centred range (from -10 to 10), a positive range (from 0 to 10), or a negative range (from -10 to 0). The mean cost per unit time as a function of γ corresponding to these three possible ranges is plotted in Fig. A.10. For a centred covariate range, we obtained the (plain) curve shown in Fig. 2, with a maximum cost for γ = 0. However, the curves differed significantly for the non-centred ranges of covariates. For a positive covariate (dashed curve), the cost increased with γ, and conversely for a negative covariate (dotted curve). This is easily understandable because, for a positive covariate, a Appendix B. Computation of threshold in case of minimal corrective repairs and dynamic covariate

This appendix presents some technical details of the computation of the threshold s λ on the conditional failure intensity in the case of minimal corrective repairs in the presence of a dynamic covariate.

As discussed in Sections 2.3 and 3, the computation of the threshold s λ requires the estimation of the function B = tφ(t) -Φ(t). Because the repairs are minimal, the estimation procedure for B (and therefore, B -1 ) can be simplified with respect to that proposed by [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF] (see Algorithm 1). The principle of this simplification is to estimate φ(t) and Φ(t) directly from the failure intensity and covariate path using the fact that Φ(t) = E[Λ t ], and not, as in [START_REF] Gilardoni | Dynamics of an optimal maintenance policy for imperfect repair models[END_REF], from M failure process simulations (using the fact that Φ(t) = E[N t ]). To this end, we simulate M covariate paths, denoted as {X (j) t } t≥0,1≤j≤M , for which we use λ (j) t and Λ (j) t to denote the corresponding intensity and cumulative intensity, respectively. In the case of minimal repair, we can estimate φ(t) and Φ(t) directly from simulations of the covariate paths:

and

where we use the notation

Finally, B(t) can be estimated as B(t) := t φ(t)-Φ(t) and we can then compute ŝλ = φ( B-1 (C P rev /C Corr )).

To achieve this, we must first determine in which time interval (denoted as l (t) in (B.2)) the inverse of the function B must be computed. As B increases, we may compute B(ldt cov ) iteratively at the end of each time interval from l ≥ 1 and stop as soon as this quantity exceeds C P rev /C Corr :