E Mosayebi Omshi

A Grall

S Shemehsavar

A dynamic auto-adaptive predictive maintenance policy for degradation with unknown parameters

Keywords: Maintenance, Condition-based, Remaining useful life, Increasing stochastic process, Bayesian update

With the development of monitoring equipment, research on condition-based maintenance (CBM) is rapidly growing. CBM optimization aims to find an optimal CBM policy which minimizes the average cost of the system over a specified duration of time. This paper proposes a dynamic auto-adaptive predictive maintenance policy for single-unit systems whose gradual deterioration is governed by an increasing stochastic process. The parameters of the degradation process are assumed to be unknown and Bayes' theorem is used to update the prior information. The time interval between two successive inspections is scheduled based on the remaining useful life (RUL) of the system and is updated along with the degradation parameters. A procedure is proposed to dynamically adapt the maintenance decision variables accordingly. Finally, different possible maintenance policies are considered and compared to illustrate their performance.

Introduction

Maintenance decision making is of prime importance for improving the global performance of industrial systems and structures during their useful life. To date, many maintenance and replacement models, from simple age-based to more sophisticated condition-based, have been proposed. Wang [START_REF] Wang | A survey of maintenance policies of deteriorating systems[END_REF] made a clear classification of different possible strategies. The tendency towards condition-based maintenance (CBM) policies is mainly due to a great improvement in technology which enables accurate online measurements of degradation levels and leads to benefits over traditional time-based maintenance [START_REF] De Jonge | The influence of practical factors on the benefits of condition-based maintenance over time-based maintenance[END_REF].

A recent review of existing maintenance policies can be found in the paper of Alaswad and Xiang [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] which also classifies CBM policies depending on the deterioration model. The major part of their paper discussed maintenance strategies in which the degradation evolution is described by stochastic processes; the Wiener, the gamma and the inverse Gaussian (IG) processes received special mention. Elwany et al. [START_REF] Elwany | Structured replacement policies for components with complex degradation processes and dedicated sensors[END_REF], Guo et al. [START_REF] Guo | A maintenance optimization model for missionoriented systems based on Wiener degradation[END_REF], and Zhang et al. [START_REF] Zhang | Degradation-based maintenance decision using stochastic filtering for systems under imperfect maintenance[END_REF] investigated examples of CBM policies based on the Wiener process. A vast amount of literature is devoted to the use of the gamma process in maintenance modeling; see, e.g., Grall et al. [START_REF] Grall | Continuous-time predictivemaintenance scheduling for a deteriorating system[END_REF], Dieulle et al. [START_REF] Dieulle | Sequential condition-based maintenance scheduling for a deteriorating system[END_REF],

Tan et al. [START_REF] Tan | Condition-based maintenance policy for gamma deteriorating systems[END_REF], Mercier et al. [START_REF] Mercier | On the modelling of imperfect repairs for a continuously monitored gamma wear process through age reduction[END_REF], and van Noortwijk [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]. The IG process has been considered more recently in degradation modeling. Chen et al. [START_REF] Chen | Condition-based maintenance using the inverse Gaussian degradation model[END_REF] investigated the optimal CBM policy with periodic inspections when the system degradation follows an IG process with a random-drift model, while Li et al. [START_REF] Li | An optimal condition based maintenance strategy using inverse-Gaussian degradation process[END_REF] obtained the optimal CBM strategy by maximizing the availability under cost constraints.

A common, but usually false assumption in maintenance literature is that the parameters of the model are known. In fact the model parameters are usually unknown and must be estimated, see, e.g., [START_REF] Chen | Condition-based maintenance using the inverse Gaussian degradation model[END_REF][START_REF] Kallen | Optimal maintenance decisions under imperfect inspection[END_REF][START_REF] Zhang | Degradation-based maintenance decision using stochastic filtering for systems under imperfect maintenance[END_REF]. De Jonge et al. [START_REF] De Jonge | Cost benefits of postponing timebased maintenance under lifetime distribution uncertainty[END_REF] considered time-based preventive maintenance under uncertainty in lifetime distribution and analyzed the effect of postponing maintenance actions on cost. The case when the degradation parameters change after a change in the system degradation rate and the new parameters are unknown was considered by Fouladirad and Grall [START_REF] Fouladirad | On-line change detection and condition-based maintenance for systems with unknown deterioration parameters[END_REF]. Kallen and van Noortwijk [START_REF] Kallen | Optimal maintenance decisions under imperfect inspection[END_REF] employed the Bayes method to estimate the unknown degradation parameters.

Articles on iterative Bayesian updates of degradation parameters can be found in Gebraeel et al. [START_REF] Gebraeel | Residual-life distributions from component degradation signals: A Bayesian approach[END_REF], Xu and Wang [START_REF] Xu | Rul estimation using an adaptive inverse gaussian model[END_REF], and Bian and Gebraeel [START_REF] Bian | Computing and updating the first-passage time distribution for randomly evolving degradation signals[END_REF]. Examples of adaptive preventive maintenance policies using Weibull distribution with the help of Bayes' theorem can be found in [START_REF] Juang | A Bayesian method on adaptive preventive maintenance problem[END_REF][START_REF] Sheu | A Bayesian approach to an adaptive preventive maintenance model[END_REF]. In the framework of dynamic condition-based maintenance, Flage et al. [START_REF] Flage | Safety constraints applied to an adaptive Bayesian condition-based maintenance optimization model[END_REF] gave a general description of the case of unknown degradation parameters with Bayesian updates and optimization of decision variables after inspections. The authors introduced and emphasize on safety constraints in decision variable optimization throughout the paper. In their study, they employed a discrete stochastic process for modeling degradation which leads to i-fold integrals for calculating the maintenance cost.

In the same framework as in [START_REF] Flage | Safety constraints applied to an adaptive Bayesian condition-based maintenance optimization model[END_REF] and to avoid overwhelming calculations, we propose a heuristic way to assess the maintenance cost. The present work considers a system subject to degradation when the degradation evolution follows a non-decreasing continuous-state stochastic process. The main contribution of the present work is to propose an aperiodic condition-based maintenance policy when the parameters of the model are unknown. All the information about these parameters (expert judgments or previous tests) may be modeled by a prior distribution. The RUL function allows us to assume such aperiodic inspection planning beside to having a high reliability (safety) level. The special case of IG process is developed because of analytical tractability in a Bayesian framework. We show how to employ Bayes' theorem to obtain successive updates for the posterior distribution of the parameters, thus increasing knowledge of the model parameters with inspection results. We also study the effect of these updates on the RUL function, used in inspection planning, as well as on the behavior of the cost function.

The remainder of the paper is organized as follows. Section 2 provides a basic description of the system and states the main assumptions of the maintenance policy. A static aperiodic maintenance policy is studied in Section 3 for the cases of both known and unknown model parameters. Section 4 describes the way the procedure of Bayesian updating is considered in adapting current information and making the maintenance decision rule dynamic and auto-adaptive in case of unknown degradation parameters. Section 5 gives more insight related to an implementation with the IG degradation process. In order to demonstrate the validity of the dynamic auto-adaptive predictive maintenance policy, different maintenance policies are compared in Section 6. Finally, Section 7 contains a conclusion.

System and Maintenance Decision Rule: Descriptions and Main Assumptions

Consider a single component system subject to wear or deterioration. Let X t denote the degradation state of the system at time t. In the absence of repair or replacement actions, the evolution of the system deterioration is assumed to be strictly increasing. This property and the time dependency structure of X t is a major motivation for employing a monotone increasing stochastic process to model the degradation evolution. We shall assume that the initial degradation state of the system is zero; i.e. X 0 = 0.

It is assumed that continuous monitoring is not possible for technological or economic reasons; the system condition can only be observed by inspection. Moreover, it is also assumed that inspec-tions are instantaneous, perfect and non-destructive. The system fails when its degradation level reaches or crosses a pre-specified threshold level L, chosen according to the system characteristics.

A failure is not self-announcing and can only be detected by inspection.

Stochastic Degradation Process

We assume that the stochastic degradation process X t can be modeled by a monotone increasing stochastic process such that 1. X 0 = 0 with probability one; 2. Non-overlapping increments of X t are statistically independent;

3. The process X t has stationary increments i.e. the distribution of X(t+s)-X(s) depends only on the length t of the interval but not on s. We refer to the pdf (resp. cdf) of the distribution of X t by f Θ(t) (x) (resp. F Θ(t) (x)). The function Θ is the vector of distribution parameters.

For the sake of simplicity, when t is omitted Θ will denote the vector of parameters which characterizes the function Θ(t).

As mentioned above, once degradation exceeds the threshold L, failure occurs and corrective action is necessary. Hence the lifetime of a product is simply the time of first passage of X t across the threshold L. At a given time t, this lifetime is fully determined by the distribution of the time until the degradation signal reaches the failure threshold, given the knowledge of the condition of the system at time t. Hence, given the parameters, we define the RUL R t,xt of a system at time t with observed degradation x t as the random variable:

R t,xt = inf{r > 0 : X t+r ≥ L|X t = x t }.
The cdf of R t,xt given the known parameters, will be denoted by F R (r|Θ, X t = x t ). Since the stochastic process is increasing, it can be expressed as:

F R (r|Θ, X t = x t ) = P (R t,xt ≤ r|X t = x t ) = P (X t+r ≥ L|X t = x t ) = P (X t+r -X t ≥ L -x t ) = F Θ(r) (L -x t ).
(1)

Predictive Maintenance Decision Rule

Let {T n } n∈N be the aperiodic sequence of inspection times (T 0 = 0). At each inspection, one must take the required maintenance decision according to the condition of the system at that time. To avoid the occurrence of system failure, a preventive threshold M < L is chosen such that preventive action is taken when the degradation level crosses this threshold. The value of M is a decision variable to be optimized. We assume that the maintenance actions are performed in a negligible time and T - n refers to the time just before the maintenance. The possible scenarios which can arise are:

• If X T -
n ≥ L, the system has failed and is correctively replaced.

• If M ≤ X T - n < L, the system has not yet failed but it has deteriorated to the extent that it is considered as worn out. In this case a preventive action is taken.

• If X T - n < M , the system is still properly functioning and there is no need for replacement. The system is left as it is.

Both preventive and corrective replacements are perfect and reset the system to "as good as new".

Hence after the inspection we have:

X Tn =    0 if X T - n ≥ M, X T - n if X T - n < M.
In all cases, the inspection scheduling is carried out based on the RUL which is defined as the remaining duration for which the system will work before it fails. The main idea of an RUL based inspection plan is that the next inspection time is chosen so that the probability of failure before the inspection remains lower than a value p (0 < p < 1). Examples of such inspection plannings can be found in [START_REF] Khoury | A comparison of RUL-based and deteriorationbased maintenance policies for gradually deteriorating systems[END_REF] and more recently in [START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF]. For this reason, let us define the p-quantile of the RUL distribution:

τ p (X Tn ) = { t : Pr(X Tn+ t ≥ L|X Tn ) = p} = { t : F R ( t|Θ( t); X Tn ) = p}. (2)
Then the next inspection time is given by

T n+1 = T n + τ p (X Tn ).
This inspection plan provides a reliability (safety) level equal to (1 -p) in which p is a decision variable to be jointly optimized with M . We remark that periodic inspections can be considered by replacing τ p (X Tn ) by a constant throughout the paper. In that case, this constant becomes a decision variable instead of p.

Cost Function

The inspections are planned sequentially and each of them incurs a cost C i . At each inspection, a preventive or corrective action is performed if necessary, with costs C p and C c respectively. Clearly, C c > C p . Moreover, since failure can only be detected through inspection, there is a system downtime after failure and an additional cost at a rate of C d is incurred from the failure time until the next replacement time. Hence the cumulative maintenance cost is:

C(t) = C i N i (t) + C p N p (t) + C c N c (t) + C d d(t),
where N i (t), N p (t), and N c (t) are respectively the number of inspections, the number of preventive replacements, and the number of corrective replacements in [0, t]. Furthermore, d(t) is the total time passed in a failed state in [0, t). In other words:

N i (t) = n∈N I {t≥Tn} , N p (t) = n∈N I {t≥Tn} I {M ≤X T - n <L} , N c (t) = n∈N I {t≥Tn} I {X T - n ≥L} , d(t) = t 0 I {Xs≥L} ds.
The ratio of the expected cost on one regenerative cycle over the expected length of the regenerative cycle will be considered in the following to determine the decision variables. Let us define it as:

EC = E[C(S 1 )] E[S 1 ] ,
where S 1 is the length of one renewal cycle i.e. the time between two successive replacements.

From classical renewal theory, it is well known that if the degradation process has the regenerative property, EC is equal to the long-run average cost per time unit EC ∞ for almost any realization of the process, i.e.,

EC = E[C(S 1 )] E[S 1 ] = lim t→∞ E[C(t)] t .
Expectations over a renewal cycle are tricky in case of aperiodic inspections. Hence, EC is obtained from EC ∞ by considering Markov renewal properties as explained in next section.

Static Aperiodic Maintenance Policy: System State and Long Run Maintenance

Cost

In this section, we characterize the evolution of the maintained system at steady state and derive an expression for the long run expected maintenance cost. The considered maintenance policy is based on the maintenance decision rule given in Section 2.2. It is qualified as "static" which means that the decision variables p and M do not change over time. They are determined in order to minimize the cost ratio EC. The following subsection is related to the classical case with known degradation parameters. The extension to unknown degradation parameters is given in the second subsection.

The Case of Known Degradation Parameters

In this subsection, we discuss the case in which the vector of parameters Θ of the degradation process is known. Let {Y n = X Tn } n∈N be the discrete-time random process describing the system state at each inspection time. Grall et al. [START_REF] Grall | Continuous-time predictivemaintenance scheduling for a deteriorating system[END_REF] derived the properties of the process {X t } t≥0 and the embedded chain {Y n } n∈N , when {X t } t≥0 is a gamma process. These properties are extended and listed for any monotone increasing stochastic process as follows. We consider the time between inspections to be based on the p-quantile of the system's remaining useful life, i.e., on τ p (•) (see ( 2)). It is worthwhile mentioning that the case of periodic inspection can be obtained by replacing the function τ p by a constant. All the calculations remain the same.

• The process {X t } t≥0 is a regenerative process with regeneration times being the dates of replacement. Also, it is a semi-regenerative process with semi-regeneration times being the inspection times.

• The discrete time process {Y n } n∈N is a Markov chain with continuous state space [0, M ) and if δ 0 denotes the Dirac delta function, the transition probability density function (conditioned on the current state) is:

Pr(dy|x, Θ) = F Θ(τp(x)) (M -x)δ 0 (dy) + f Θ(τp(x)) (y -x)I {x≤y<M } dy, (3)
where

F Θ(τp(x)) (y) = +∞ y f Θ(τp(x)) (u)du.
The probability function Pr(dy|x, Θ) describes the probability of transition from state x to state y. The mass part is related to the degradation level y = 0 which is obtained after maintenance. The density part describes the degradation in the absence of maintenance.

• The Markov chain {Y n } n∈N is Harris ergodic with a unique stationary probability distribution

π Θ : π Θ (dx) = a Θ δ 0 (dx) + (1 -a Θ )b Θ (x)dx (4)
with

a Θ = 1 1 + M 0 B Θ (x)dx , b Θ (y) = a Θ 1 -a Θ B Θ (y),
where

B Θ (y) = f Θ(τp(0)) (y) + y 0 B Θ (x)f Θ(τp(x)) (y -x)dx.
• The following result follows from the fact that {X t } t≥0 is a semi-regenerative process and the unique stationary probability distribution π Θ exists for the embedded Markov chain {Y n } n∈N .

We have:

EC = E[C(S 1 )] E[S 1 ] = EC ∞ = E πΘ [C(T 1 )] E πΘ [T 1 ] ,
where T 1 is the first inspection time. This means that an expression for the maintenance cost EC, as well as EC ∞ , can be obtained from the analysis of the system behavior on a semiregenerative cycle (i.e. between two inspections) instead of a regenerative cycle (i.e. between two replacements).

Therefore, we have:

EC = C i E πΘ [N i (T 1 )] E πΘ [T 1 ] + C p E πΘ [N p (T 1 )] E πΘ [T 1 ] + C c E πΘ [N c (T 1 )] E πΘ [T 1 ] + C d E πΘ [d(T 1 )] E πΘ [T 1 ] . (5)
Obviously,

E πΘ [N i (T 1 )] = 1. The other quantities, E πΘ [N p (T 1 )], E πΘ [N c (T 1 )], E πΘ [d(T 1 )],
and E πΘ [T 1 ], are computed by integration with respect to π Θ as follows:

E πΘ [N p (T 1 )] = P πΘ (M ≤ X T - 1 < L) = M 0 F Θ(τp(x)) (M -x) -F Θ(τp(x)) (L -x) π Θ (dx), E πΘ [N c (T 1 )] = P πΘ (X T - 1 ≥ L) = M 0 F Θ(τp(x)) (L -x)π Θ (dx), E πΘ [d(T 1 )] = M 0 τp(x) 0 F Θ(τp(x)) (L -x)ds π Θ (dx), E πΘ [T 1 ] = M 0 τ p (x)π Θ (dx).

The Case of Unknown Degradation Parameters

In practice, the vector of parameters Θ of the degradation process is not known. In this paper, we

propose to use the available information about the parameters as a prior distribution, then update this information using Bayesian methods following the inspection planning which is given by the decision rule. Let us first define the prior distribution of Θ which depends on hyperparameters Ξ.

It is denoted by g Ξ hereafter.

The distribution of the residual life of the system is obtained from g Ξ and employed to assess the next inspection time. Knowing the prior distribution of Θ, the cdf of remaining useful life can be obtained by:

F R (r|X t = x t ) = F R|Θ (r|Θ, X t = x t )g Ξ (Θ)dΘ. ( 6
)
For unknown parameters and given the prior distribution, the four properties of the process {X t } t≥0 and of the chain {Y n = X Tn } n∈N mentioned in Section 3.1 can be restated as follows.

• The process {X t } t≥0 is a regenerative and semi-regenerative process with regeneration times and semi-regeneration times, S 1 and T 1 , respectively.

• The process {Y n } n∈N is a Markov chain that takes values in [0, M ) and has transition probability density function (conditioned on the current state)

Pr(dy|x) = F * τp(x) (M -x)δ 0 (dy) + f * τp(x) (y -x)I {x≤y<M } dy, (7)
where

f * τp(x) (z) = f Θ(τp(x)) (z)dΘ, (8)
and

F * τp(x) (z) = ∞ z f * τp(x) (u)du.
• Similarly, the unique stationary probability distribution of {Y n } n∈N is given by:

π(dx) = aδ 0 (dx) + (1 -a)b(x)dx, (9)
with

a = 1 1 + M 0 B(x)dx , b(y) = a 1 -a B(y),
where

B(y) = f * τp(x) (y) + y 0 B(x)f * τp(x) (y -x)dx.
• As for the case with known degradation parameters, we find:

EC = EC ∞ = E π [C(T 1 )] E π [T 1 ] . ( 10
)
The expectations for the different costs are computed as follows:

E π (N p (T 1 )) = P π (M ≤ X T - 1 < L) = M 0 F * τp(x) (M -x) -F * τp(x) (L -x) π(dx), E π (N c (T 1 )) = P π (X T - 1 ≥ L) = M 0 F * τp(x) (L -x)π(dx), E π (d(T 1 )) = M 0 τp(x) 0 F * s (L -x)ds π(dx), E π (T 1 ) = M 0 τ p (x)π(dx).

Optimization of The Maintenance Policy

The choice of the maintenance threshold value and inspection dates will obviously influence the economic performance of the maintenance policy. Thus, the objective is to obtain values of the decision variables p and M such that the maintenance cost EC reaches its minimum value. For the static aperiodic maintenance cases described in Sections 3.1 and 3.2 respectively, the optimal decision variables, p * and M * , are such that:

EC(p * , M * ) ≤ EC(p, M ), ∀(p, M ) ∈ [0, 1] × [0, L],
where EC(p, M ) is as in ( 5) or [START_REF] Grall | Continuous-time predictivemaintenance scheduling for a deteriorating system[END_REF], respectively. Numerical illustrations of this procedure are given in Section 6.

From Static to Adaptive Predictive Aperiodic Maintenance Policy

The previous section gives probabilistic characteristics of the system evolution under a given aperiodic maintenance policy. It allows us to derive a cost rate EC and hence optimize the decision variables p and M , which are presumed to be constant over time. In that case, the cost EC is equal to the long run expected maintenance cost per time unit. Section 4.1 describes a procedure which updates the estimates of the model parameters using Bayesian methods and adjusts the maintenance decision rule accordingly. In that case, the cost EC which is obtained at a given cycle is related to current level of information and it changes with the change of model parameters. As a consequence EC is no longer equal to EC ∞ . We heuristically propose to use the cost ratio EC as a criterion to determine the decision variables according to the current level of knowledge on the degradation process. Three alternative maintenance policies are introduced and their properties are derived in Section 4.2.

Degradation Model Update and Adaptive Maintenance Decision Rule

This Bayesian procedure allows us to update current information about Θ with a new observation. Updating can take place at each inspection or at the end of each cycle with the corresponding number of inspections. Bayesian updates can be used directly to modify the inspection scheduling function. The algorithm is as follows:

• Update after each inspection. In this case, only one new observation is available. For a given p and M we have:

Step 1: Initialize the hyperparameters Ξ (0) . Furthermore, set X 0 = 0, T 0 = 0, and i = 0.

Step 2: Using (6), evaluate t = τ p (X Ti ) according to the RUL distribution.

Step 3: Monitor the system at T i+1 = T i + t.

Step 4: Using Bayes' theorem, obtain the joint posterior distribution of Θ with the information available (i.e., X Ti and X Ti+1 ) to obtain the new hyperparameters Ξ (i+1) .

Step 5: If X Ti+1 ≥ M , i.e. the system has failed or is too deteriorated, set X Ti+1 = 0 (corrective or preventive replacement, whichever is needed).

Step 6: Set i = i + 1 and then return to Step 2.

• Update after each cycle. The only difference with the update after each inspection rule is that computation of the posterior distribution of Θ is postponed until the end of each cycle and the hyperparameters remain unchanged during a specific cycle. Here, the number of the new observations available at each update is equal to the number of inspections in that cycle.

Static and Dynamic Auto-adaptive Maintenance Policy

The maintenance decision rule described in Section 2.2 depends on two decision variables which are the value of p for the p-quantile of the RUL and the preventive threshold M . As stated in the previous section, the distribution of model parameters and hence the RUL distribution and the expression of the time to the next inspection can be updated without changing the decision variables. In this section, we consider options for maintenance policies as the values of the decision variables or the distribution of the parameters are modified.

Depending on decision variables, policies can be classified in two groups:

Static strategy: The decision variables are fixed over the time horizon in which the policy is implemented.

This strategy can be implemented when the degradation parameters are known or unknown, based on the results in Section 3. For known degradation parameters the expression of τ p is obtained according to the RUL function given in [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF]. For unknown degradation parameters the p-quantile of the RUL distribution, given in ( 6), can remain constant based on prior information or be updated as explained in Section 4.1. The corresponding maintenance policies are respectively called "static non-adaptive" or "static auto-adaptive".

Dynamic strategy:

The decision variables can be modified "on-line" and be updated from information which is collected at inspection times.

The resulting policy, which is implemented in the case of unknown degradation parameters, is called a "dynamic auto-adaptive" policy.

To summarize, three different aperiodic predictive maintenance policies can be implemented: In a general case, the derivation of optimal decision variables p and M with degradation parameters updates is a difficult problem. Consequently a generic process is proposed which is based on minimization of the maintenance cost EC given the current available information. This cost can be obtained based on analytical calculations as described in Section 3, The specific proposed rules are as follows:

• For SNAP policy implementation, the fixed decision variables are directly obtained by minimizing EC as derived respectively in Section 3.1 for known degradation parameters and in Section 3.2 for unknown parameters. In the second case the prior distribution of parameters Θ is with initial values of the hyperparameters. In that configuration, EC is equal to the long-run cost rate EC ∞ .

• For implementation of the DAPAM policy a heuristic suboptimal procedure is proposed in which the decision variables are optimized for each cycle according to the available knowledge at the beginning of the cycle. The resulting values are used for one renewal cycle, i.e. until the next update. The optimization is based on the cost rate EC as given in Section 3.2 and is conducted using the most recent updated values of the hyperparameters.

Concerning the SAAP policy, the decision variables are fixed. They should be obtained at the initial time, before any update. This policy is considered only for numerical reference in Section 6.

Monte Carlo simulation is performed with degradation levels obtained from real values of parameters to estimate the maintenance cost rate.

In the next sections, the precise implementation is given specifically for the inverse Gaussian process. Then the policies are analyzed and compared numerically using the mean estimate of EC as a criterion in order to assess the updating processes and the heuristic choice of decision variables.

This estimate is obtained on a finite time horizon.

Case of the Inverse Gaussian process

As indicated in the introduction, different stochastic processes can be considered, for example the gamma process and the IG process. Both have the same meaningful physical interpretation for modeling degradation as a limit of a compound Poisson process [START_REF] Ye | The inverse Gaussian process as a degradation model[END_REF]. In addition these two processes have a high capacity to incorporate random effects and covariates to capture heterogeneities in degradation data. Albeit, the IG process is more appropriate for incorporating prior information, especially conjugate information, into Bayesian inference due to its relationship with the Gaussian distribution family. Explicit forms of important related characteristics such as reliability function, quantiles, and the distribution function of the remaining useful lifetime are also easily obtained.

Therefore, the IG process is specifically considered in this section for the implementation of the proposed adaptive policies.

For this process, each increment follows an IG distribution, that is, X t -X s ∼ IG(δ -1 .(ts), λ.(t -s) 2 ) for all t > s. The related pdf and cdf are defined for x > 0 as follows:

f Θ(t) (x) = f δ -1 t,λt 2 (x) = λt 2 2πx 3 exp{- λ 2x (δx -t) 2 }, F Θ(t) (x) = F δ -1 t,λt 2 (x) = Φ √ λ(t -δx) √ x -exp(2tλδ)Φ - √ λ(t + δx) √ x
where δ, λ > 0 and Φ is the cdf of the standard normal distribution. Then, the mean and variance of X t are t/δ and t/(δ 3 λ), respectively. From (1), the cdf of R t,xt given the known parameters, can be expressed as:

F R|δ,λ (r|δ, λ, X t = x t ) = Φ √ λ(r -δ(L -x t )) √ L -x t -exp(2rλδ)Φ - √ λ(r + δ(L -x t )) √ L -x t .
With these considerations, the functions in (3), ( 4) and ( 5) can easily be obtained.

In the case of unknown parameters and for Bayesian inference, conjugate priors are employed to ease calculations. Hence, suppose λ has the gamma density function

f (λ) = λ α-1 Γ(α)β α exp{-λ/β}, (11)
and let δ have the conditional normal density function with mean ξ and variance σ 2 /λ:

f (δ|λ) = λ 2πσ 2 exp{- λ(δ -ξ) 2 2σ 2 }. (12)
Then the joint prior distribution of (δ, λ) is given by f (δ, λ) = f (δ|λ)f (λ). With this definition, some characteristics of δ and λ are

E(λ) = αβ, Var(λ) = αβ 2 , E(δ) = ξ, and Var(δ) = σ 2 (α -1)β .
Clearly, the slope of an irreversible degradation process must be positive. Hence δ can only take positive real values. To avoid the possibility of obtaining negative degradation slopes, we suppose that P (δ ≤ 0) is negligible (Lu and Meeker [START_REF] Lu | Using degradation measures to estimate a time-to-failure distribution[END_REF], and Peng [START_REF] Peng | Inverse Gaussian processes with random effects and explanatory variables for degradation data[END_REF]). With a calculation similar to that in Peng [START_REF] Peng | Inverse Gaussian processes with random effects and explanatory variables for degradation data[END_REF], ( 6) and ( 8) are as follows:

F R (r|X t = x t ) = 1 - β 2π Γ(α + 1/2)r Γ(α) L-xt 0 z -3/2 (σ 2 z + 1) -1/2 1 + β(ξz -r) 2 2z(σ 2 z + 1) -(α+1/2)
dz, and

f * τp(x) (z) = β 2π Γ(α + 1/2)τ p (x) Γ(α) z -3/2 (σ 2 z + 1) -1/2 1 + β(ξz -τ p (x)) 2 2z(σ 2 z + 1) -(α+1/2)
.

Then, functions in ( 7), ( 9) and ( 10) can be obtained. The posterior functions of (δ, λ) needed in Section 4.1 are also given in the Appendix.

Numerical Study

A simulation study is conducted in the case of an IG process for better insight through numerical illustration. The aim of the numerical experiments is twofold.

i) To assess the performance of the proposed auto adaptive maintenance policies and compare them with non-adaptive policies.

ii) To study the behavior of the decision variables and model parameters and highlight the effect of the Bayesian updates.

Simulation Framework

Samples of degradation increments are generated from an IG process with fixed parameters δ real and λ real . Three values of the mean degradation rate (0.6, 1 and 2) and three values of the variance (0.25, 1 and 4) are considered. The nine resulting configurations are summarized in Table 1. The failure level is set at L = 9 for all samples. In the case of unknown parameters a prior distribution is assumed as described in Section 5. The initial values considered for hyperparameters are as follows:

α = 1.5, β = 2 3 , ξ = 1, and σ = 1 √ 3 . (13)
Figure 1 displays the graph of the corresponding joint pdf of the degradation parameters. The mean and variance of the degradation parameters for this first prior are respectively increases with σ (respectively α.β 2 ). A second example of a prior distribution with more uncertainty about degradation parameters will be considered in Section 6.2.3 so as to analyse sensitivity to choice of prior. The value of σ will be significantly increased with the second prior.

E[λ] = E[δ] = 1, Var[δ] = 1,

Comparison of Maintenance policies

In order to assess and compare the different maintenance policies, hence the relevance and responsiveness of the dynamic auto-adaptive process, the behavior of the system is considered for a limited time horizon (0 to T end = 200). The criterion for comparison is the maintenance cost estimated on the chosen time horizon. Classical Monte Carlo simulation is considered for estimation: the maintenance decision rules are applied on simulated degradation paths. For the DAPAM policy, the optimization for decision variables is conducted by a classical gradient-based optimization algorithm and requires the most computational effort. With R software on an Intel core i5 processor the time required to calculate a cost value (stationary law + cost assessment) is about 30 seconds. Around 30 cost evaluations may be necessary for convergence of one optimization.

Static non-adaptive (SNAP) maintenance policy

Let consider first the case of static maintenance policy without Bayesian update. The two options with known and unknown degradation parameters are investigated hereafter. In both cases the decision variables are constant and optimized from information available at time t = 0. .

• Known degradation parameters

The maintenance decision rule is derived and tuned assuming that the degradation follows an IG process with known parameters. The values p * and M * of the decision variables are those minimizing EC in [START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF]. These values are given in the last three rows of Table 2 for each configuration. Rows 2 and 3 correspond to the values of the maintenance cost estimated on a finite time horizon EC and its standard deviation. The estimation is based on simulation of 45 samples of maintained system degradation paths over a finite time horizon.

It can be seen in Table 2 that the values of theoretical and estimated maintenance cost rates, i.e. EC and EC, are close to each other. Moreover, the standard deviations of the estimations are very small. In other words, the decision variables optimized with respect to EC are relevant with estimation on a finite time horizon.

• Unknown degradation parameters

Here, an IG process is assumed with the proposed prior distribution of hyperparameters given in [START_REF] De Jonge | Cost benefits of postponing timebased maintenance under lifetime distribution uncertainty[END_REF]. The initial values of the hyperparameters are not updated. Analytical calculations from Section 3.2 are used to obtain and minimize EC. The optimal decision variables and the corre- The values of p and M given in ( 14) can be used as decision variables to apply the SNAP policy to the configurations described in Table 1. The related estimated maintenance costs on the finite time horizon are given in Table 3. As expected, they are higher than for the case of known parameters.

However, the policy seems to be robust. Here, we assume that the degradation parameters are unknown and consider policies involving Bayesian updates. First, the case with constant decision variables (the SAAP maintenance policy) is considered, and then the case where the decision variables are automatically adapted to the update of the hyperparameters (the DAPAM policy) is considered.

• Static auto-adaptive (SAAP) maintenance policy

The hyperparameters are updated at each inspection time. The expression for the p-quantile of the RUL is modified accordingly but the values of p and M remain constant over the whole time horizon. The maintenance cost EC depends on the characteristics of the Bayesian update sequence up to T end and these are related to the real degradation parameters. Hence the optimization of the decision variables requires knowledge of the real degradation parameters. These parameters are not available in reality but numerical experiments are achievable. Table 4 gives the minimum estimated values of the cost rate over a finite horizon for the SAAP policy with updates at each inspection time. A comparison of the values in Table 4 with the values in the last three rows in Table 2 reveals the efficiency of the Bayesian update in the maintenance decision rule. The minimum values of EC are very close to those obtained for known decision variables, which means that the influence of the prior weakens over time, and that the Bayesian update is efficient from the maintenance point of view.

M * 7. 2 6.4 6.2 7.4 6.4 6.4 7.6 7 6.6
• Dynamic auto-adaptive (DAPAM) maintenance policy

The DAPAM policy allows decision variables to be optimized sequentially. The optimization can be performed at any time according to the available information described by the hyperparameters at that time.

The estimation of the maintenance cost on successive renewal cycles requires Monte Carlo simulation of degradation paths followed by optimization of p and M . The optimization involves numerical calculations for the stationary law and long-run cost assessment as described in Section 3.2. The expected maintenance cost on a finite time horizon EC is estimated after the last update when time T end has been reached. The results are given in Table 5. We see that the cost values are close to those of the SNAP policy with known parameters (Table 2) for each configuration. The relative deviations between the costs of the two policies are lower for high degradation rates (e.g δ = 0.5) than for low degradation rates (e.g. δ = 1.667). This is because the number of renewal cycles and hence the number of updates from 0 to T end increases with the mean deterioration speed.

The smallest relative difference of costs is obtained for δ = 1 which is the case for which the initial values of hyperparameters fit best. Overall, the results show the promise of the DAPAM policy. In the case of unknown degradation parameters, the initial knowledge about the degradation may have an influence on the economic performance of the auto-adaptive policies. In order to illustrate this, let us consider another set of initial values for hyperparameters as follows:

α = 3, β = 2 3 , ξ = 1 2 , and σ = 12. (15)
The mean and variance of the degradation parameters for this prior are respectively The graph of the pdf of (λ, δ) is displayed in Figure 2. For these values of the hyperparameters, the SNAP policy for the optimization of the long term cost rate leads to p * = 0.603, M * = 6.24, and EC * = 2.43. [START_REF] Khoury | A comparison of RUL-based and deteriorationbased maintenance policies for gradually deteriorating systems[END_REF] This cost can be compared with the cost in [START_REF] Juang | A Bayesian method on adaptive preventive maintenance problem[END_REF] which is significantly lower. Actually, the distribution of degradation parameters from [START_REF] Guo | A maintenance optimization model for missionoriented systems based on Wiener degradation[END_REF] and ( 12) is less informative with the set of hyperparameters given in [START_REF] Kallen | Optimal maintenance decisions under imperfect inspection[END_REF] than with the set of hyperparameters given in (13) (the variances of δ and λ are greater). It follows that the density function of (λ, δ) is more spread out. The numerical results corresponding to the different policies are given in Table 6 for Configurations 5 and 9. The cost standard deviations for both cases are around 0.15. The cost values for SNAP, SAAP, and DAPAM policies can be compared with those in Tables 3, 4 and 5. The SNAP policy has been implemented assuming unknown parameters (hyperparameters from ( 15)) and decision variables as in [START_REF] Khoury | A comparison of RUL-based and deteriorationbased maintenance policies for gradually deteriorating systems[END_REF]. The SAAP policy is implemented for the best case i.e. the decision variables are those minimizing EC. 6)) is very close to the cost for the case of known parameters. This shows that the first few Bayesian updates and especially the very first one, are fundamental from the maintenance point of view. At the end of the time horizon, the behavior of the system with the DAPAM policy is very close to the behavior in the case when the degradation parameters are known. We remark that from a practical point of view it can be difficult to choose initial hyperparameters in some cases. If maintenance planners want to avoid failures and start with a conservative point of view it may be interesting to consider a prior distribution with an underestimated value of δ, i.e., with overestimation of the degradation rate. The drawback to this may be a significant increase of maintenance cost for the first cycle.

Parameters of the adaptive maintenance policy

In this section, we present three typical cases which depict the evolution of the distribution of the unknown parameters during the process of Bayesian updating for systems maintained with the DAPAM policy. With the initial values of hyperparameters as in [START_REF] De Jonge | Cost benefits of postponing timebased maintenance under lifetime distribution uncertainty[END_REF], Configurations 3 and 8 with δ = 1.667 and δ = 0.5 respectively represent examples of low and high degradation rates while Configuration 5 with initial values of hyperparameters as in [START_REF] Kallen | Optimal maintenance decisions under imperfect inspection[END_REF] illustrates the case of a poor prior. given in Table 2. We note that in the case of Configuration 5 (Figure 7), the first Bayesian update which is taken around t = 25 has a significant impact on the rate of increase of the trend. , decrease, reflecting the high certainty of the tendency. This procedure is more obvious in Figure 6. This is due to the higher number of updates which is the consequence of the higher degradation rate in configuration 8. Once again, Figure 8 shows the influence of the Bayesian updates in the convergence even when the priors have not been very well chosen.

Conclusion

Maintenance has become an increasingly important area of research and a large body of literature has been developed for it. However, almost all studies share the assumption that the parameters of the degradation process are known. In this paper, we consider a dynamic auto-adaptive predictive maintenance (DAPAM) policy for single-unit systems with unknown degradation parameters, assuming that the degradation is governed by a monotone increasing process and that the unknown parameters jointly follow a prior distribution with specified hyperparameters. The Bayes method is employed to update information about parameters over time. Simulation studies are conducted to illustrate the behavior of the policy and maintenance decision variables are obtained using a mean cost rate between two replacements as a criterion in the case of degradation governed by an IG process. Moreover, different alternative policies have been discussed and it is shown that when compared with the DAPAM policy the latter performs much better when considering its ability to get close to the optimal decision variables obtained in the case that the degradation parameters are known.

The proposed dynamic decision rule is heuristic. In case of unknown or uncertain parameters and for continuous deterioration processes, the optimal policy structure is unknown. Future research needs to be done to characterize properties of the proposed policy in terms of optimality or suboptimality. One next step would be to obtain analytical results concerning convergence rates of the Bayesian estimation and the long-run cost rate following given numbers of Bayesian updates.

Although some advantages of the DAPAM policy are its small number of decision variables and its adaptability to a wide range of situations, a possible weakness could be a slow convergence rate in parameters estimation in certain cases. Actually the considered inspection policy tends to reduce the number of inspections. More numerical analysis needs to be carried out in order to characterize disadvantageous cases. The influences of several things need to be analyzed more deeply including the initial prior distribution, the real degradation process, the inspection cost, and the number of inspections.

There are potential extensions to overcome the possible drawbacks concerning estimation when it is driven by the proposed condition-based inspection policy. An additional rule could be considered to increase the number of inspections when necessary. Frequent periodic inspections could be planned during a given period of time to get more information about degradation and accelerate parameters estimation. In this case a specific additional decision rule needs to be defined, with its additional decision variables, and compared to the proposed one. Natural extensions of the present work could also consider more complex degradation processes with the influence of environmental parameters as well as imperfect maintenance actions.

Appendix

The appendix gives a proposition which can be used to find the successive updates for parameter distributions in the case of an IG process. Suppose X t0 , X t1 , . . . , X t k are the observed degradation data at times t 0 , t 1 , . . . , t k . Then the likelihood function of δ and λ is:

L(δ, λ) ∝ λ k 2 exp{- λ 2 k i=1 ( x i δ -t i ) 2 x i },
where x i = X ti -X ti-1 and t i = t i -t i-1 . Using ( 11) and ( 12), The joint posterior distribution of (δ, λ) can be derived as follows:

Proposition: Given the observed data, Data = {X t0 , X t1 , . . . , X t k }, the joint posterior distribution of (δ, λ) is: where

f (δ, λ|Data) = f (
A = k i=1 x i + 1 σ 2 , B = k i=1 t i + ξ σ 2 , C = k i=1 ( t i ) 2 x i + ξ 2 σ 2 , D = 1 2 (C - B 2 A
).

Proof : The joint posterior of the (δ, λ) is f (δ, λ|Data) ∝ L(δ, λ)f (δ|λ)f (λ)

∝ λ k+1 2 +α-1 e -λ β exp{- λ 2 [ (δ -ξ) 2 σ 2 + k i=1 (δ x i -t i ) 2 x i ]} ∝ λ k+1 2 +α-1 e -λ β exp{- λ 2 A[δ 2 -2δ B A + C A ]} ∝ λ k+1 2 +α-1 e -λ β exp{- λ 2 A(δ - B A ) 2 - λ D } ∝ λ 1 2 e -λ 2 A(δ-B A ) 2 λ k 2 +α-1 e -λ[ 1 β + 1 D ] ,
where A, B, C, and D are as in the proposition. This implies that the posterior distributions of δ|λ and λ are respectively the normal and gamma distributions with the hyperparameters as given.

  (i) a Static Non-adaptive Aperiodic Predictive (SNAP) policy with fixed decision variables and degradation parameters or hyperparameters, (ii) a Static Auto-adaptive Aperiodic Predictive (SAAP) policy where only the parameters' distribution (hence the RUL distribution) is updated, (iii) the Dynamic Auto-adaptive Predictive Aperiodic Maintenance (DAPAM) policy where the parameters' distribution and the decision variables are updated simultaneously.

  and Var[λ] 0.67. The hyperparameter relationships have the following interpretations: as E[δ] = ξ, the mean degradation rate is increasing as ξ decreases. The variance of the degradation rate increases as λ decreases hence when the product α.β of the hyperparameters decreases. The variance of δ (respectively λ)

Figure 1 :

 1 Figure 1: The joint density function of (λ, δ) for α = 1.5, β = 2 3 , ξ = 1, and σ = 1 √ 3

  sponding minimum cost obtained by numerical methods are respectively: p * = 0.115, M * = 6.49, and EC * = 1.19.[START_REF] Juang | A Bayesian method on adaptive preventive maintenance problem[END_REF]

Figure 2 :

 2 Figure 2: The joint density function of (λ, δ) for α = 3, β = 2 3 , ξ = 1 2 , and σ = 12.

Figures 3 ,

 3 Figures 3,[START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF] and 7 display the evolution of the means and standard deviations of the decision variables as a function of time for these three cases. For example, Figure3shows the evolution of p and M for Configuration 3. We see that over time, the standard deviations of these estimators decrease. Their values tend to the values of p * and M * obtained in the case of known parameters as

Figures 4 ,

 4 Figures 4, 6 and 8 show the evolution of the hyperparameters. The joint distribution of (δ, λ) is such that E(δ) = ξ and E(λ) = αβ. The red lines on the graphs of the evolutions of αβ and ξ correspond to the values of δ real and λ real in the configuration. These figures show that even though δ and λ are random variables, the means of their distributions, i.e ξ and αβ, tend to the values δ real and λ real with the successive Bayesian updates. At the same time, the variance of λ and δ, i.e. αβ 2 and σ 2 (α -1)β , decrease, reflecting the high certainty of the tendency. This procedure is

Figure 3 : 1 √ 3 .

 313 Figure 3: The evolution of decision variables p (left) and M (right) for the DAPAM policy in Configuration 3 with α = 1.5, β = 2 3 , ξ = 1, and σ = 1 √ 3.

2 Figure 4 : 1 √ 3 .

 2413 Figure 4: The evolution of hyperparameters for the DAPAM policy in Configuration 3 with α = 1.5, β = 2 3 , ξ = 1, and σ = 1 √ 3.

Figure 5 : 1 √ 3 . 2 Figure 6 : 1 √ 3 .

 5132613 Figure 5: The evolution of decision variables p (left) and M (right) for the DAPAM policy in Configuration 8 with α = 1.5, β = 2 3 , ξ = 1, and σ = 1 √ 3.

Figure 7 :

 7 Figure 7: The evolution of decision variables p (left) and M (right) for the DAPAM policy in Configuration 5 with α = 3, β = 2 3 , ξ = 1 2 , and σ = 12.

Figure 8 :

 8 Figure 8: The evolution of hyperparameters for the DAPAM policy in Configuration 5 with α = 3, β = 2 3 , ξ = 1 2 , and σ = 12.

  δ|λ, Data)f (λ|Data), wheref (δ|λ, Data) = λ 2πσ * 2 exp{-λ(δ -ξ * ) 2 2σ * 2 }, f (λ|Data) = λ α * -1 Γ(α * )β * α * exp{-λ/β * }.The updated hyperparameters are ξ * = B/A, σ * = A -1/2 , α * = α + k/2, and β * = (1/β + 1/D) -1 ;

Table 1 :

 1 Different configurations for parameters of degradation process.

	Config. nb. δ real	λ real Variance
	1	1.667 0.864	0.25
	2	1.667 0.216	1
	3	1.667 0.054	4
	4	1	4	0.25
	5	1	1	1
	6	1	0.25	4
	7	0.5	32	0.25
	8	0.5	8	1
	9	0.5	2	4

Table 2 :

 2 The optimal decision variables and costs for the SNAP policy when the parameters are known.

	Config. nb.	1	2	3	4	5	6	7	8	9
	EC	0.386 0.472 0.527 0.605 0.739 0.827 1.141 1.323 1.540
	SD	0.047 0.074 0.108 0.047 0.087 0.129 0.057 0.088 0.141
	EC *	0.394 0.499 0.549 0.609 0.739 0.856 1.136 1.313 1.552
	p *	0.017 0.033 0.046 0.015 0.033 0.056 0.012 0.028 0.061
	M *	7.05	6.28	6.01	7.38	6.59	6.21	7.74	7.02	6.37

Table 3 :

 3 Values of EC and their standard deviations for the SNAP policy assessed based on the hyperparameters given in[START_REF] De Jonge | Cost benefits of postponing timebased maintenance under lifetime distribution uncertainty[END_REF].

	Config. nb.	1	2	3	4	5	6	7	8	9
	EC	0.462 0.530 0.559 0.718 0.791 0.869 1.365 1.450 1.579
	SD	0.098 0.130 0.152 0.105 0.122 0.154 0.129 0.136 0.164
	6.2.2. Auto-adaptive maintenance policies

Table 4 :

 4 Values of the optimal estimated maintenance cost over a finite time horizon and corresponding decision variables for the SAAP policy with update at each inspection time.

	Config. nb.	1	2	3	4	5	6	7	8	9
	EC	*	0.394 0.488 0.533 0.612 0.736 0.838 1.144 1.322 1.544
	p *		0.028 0.046 0.062 0.028 0.038 0.064 0.06 0.036 0.07

Table 5 :

 5 Values of the maintenance cost over a finite time horizon and corresponding decision variables for the DAPAM policy with update at each cycle.

	Config. nb.	1	2	3	4	5	6	7	8	9
	EC	0.375 0.504 0.582 0.614 0.753 0.823 1.190 1.318 1.580
	SD	0.035 0.086 0.135 0.043 0.095 0.121 0.038 0.068 0.142
	p at T end	0.016 0.036 0.050 0.018 0.036 0.058 0.020 0.032 0.068
	M at T end	7.13	6.17	5.95	7.18	6.46	6.06	7.24	6.82	6.23
	6.2.3. Effect of prior knowledge

Table 6 :

 6 Costs of EC for the SNAP, SAAP, and DAPAM policies in Configurations 5 and 9 with the hyperparameters given in Eq.[START_REF] Kallen | Optimal maintenance decisions under imperfect inspection[END_REF]. The bracketed values are accumulated costs from the beginning of the second renewal cycle toT end .This is very disadvantageous when using the DAPAM policy because it leads to a very late first inspection. Hence there is only one inspection in the first renewal cycle which is at t = 25. A corrective maintenance is performed at that time after a long period of unavailability. Due to the unsuitable initial prior distribution, the cost incurred during the first renewal cycle is very high and more than 12% of the time horizon elapses without any update. Without the cost of the first renewal cycle or, in other words, if we start counting after the first update, the cumulated maintenance cost rate (bracketed values in Table

	Config. nb.	EC	EC	EC	p at T end	M at T end
		(SNAP) (SAAP) (DAPAM) (DAPAM) (DAPAM)
	5	1.53	0.777	1.02	0.035	6.50
				(0.739)
	9	3.02	1.57	1.80	0.066	6.25
				(1.53)
	With initial values of hyperparameters given in (15), the RUL prediction is an over-estimate.


Jordan Gauvrit

Antoine Caucheteux

Stéphane Lecoeuche

Uncertainty Propagation of Internal Heat Gains for Building Thermal Behaviour Assessment: Influence of Spatial Distribution

In building simulation, internal heat gains correspond to heat production by human metabolism or electrical devices use. It is one of the most uncertain model inputs and could have an important impact on building simulation results. This study proposes a method to investigate the influence of the internal heat gains uncertainties by separating the uncertainty on the internal heat gains of the entire building, the uncertainty on the spatial distribution and its evolution on time. The uncertainty sources are propagated independently in a dynamic thermal simulation (DTS). The temperature of each zone at each moment is analyzed. In order to simplify this study, the most representing temperatures are selected with a method based on cumulative variances and a clustering algorithm. This approach is applied on an office building in France. The data coming from a one year monitoring period, provide information to reduce the uncertainties about the real internal heat gains. The results indicate that the effects of the internal heat gains uncertainties are time dependent. They also depend on the heating scenario of the thermal zone (heated or not-heated). At last, the temperatures are mainly influenced by the uncertainty on the internal heat gains of the entire building.

Introduction

Many studies have attempted to reduce the gap between the predicted and measured energy consumption by the calibration of DTS model in the building's monitoring conditions [START_REF] Cheng | The gap between predicted and measured energy performance of buildings: A framework for investigation[END_REF]. Occupant behavior is one of the key parameter in this exercise [START_REF] Mahdavi | User Behavior And Energy Performance In Buildings[END_REF]. Indeed, it is necessary to take into account the actions of the occupant on some building components (the settings of Heating Ventilation and Air Conditioning, the opening of windows…) and the internal heat gains due to the occupant body metabolism and the heat provided by the electrical devices. This study will focus on the configuration of internal heat gains for the calibration of DTS. Data about the internal heat gains of the entire building and its spatial distribution between each zone are needed. For instance, electrical meters may give information about the heat provided by the electrical devices to the entire building. Moreover, a building occupancy schedule can be used to assess the heat gain by human metabolism. These sources of information do not allow to assess precisely the internal heat gains, but various constraints as the cost of the measurement equipment or the respect of privacy limit the availability of more accurate data [START_REF] Yan | IEA EBC Annex 66: Definition and simulation of occupant behavior in buildings[END_REF][START_REF] Zhang | Rethinking the role of occupant behavior in building energy performance: A review[END_REF]. As a results, the internal heat gains of the entire building, the spatial repartition between the thermal zones and the time dependency of spatial repartition still uncertain model parameters. Therefore, the simulated results will be affected by this imprecise setting. For this reason, uncertainty propagation (UP) techniques have been applied on DTS model to assess the influence of parameters uncertainties on the model outputs. [START_REF] Ligier | Energy Performance Contracting Methodology Based upon Simulation and Measurement[END_REF] or [START_REF] Ren | The Influence of the Variation in Occupancy Patterns on Domestic Energy Simulation Prediction: A Case Study in Shanghai[END_REF] have assessed the effects related to the uncertainty on the internal heat gains of the entire building. In this work, the uncertainty due to the spatial repartition and its time dependency are also taken into account. The uncertainty sources have been propagated with a Monte-Carlo procedure of simulations [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF]. [START_REF] Bipm | Evaluation of measurement data -Supplement 1 to the "Guide to the expression of uncertainty in measurement" -Propagation of distributions using a Monte Carlo method[END_REF] recommends to choose the number of simulations such as the confidence level of the model output interval reaches a desired value. For a 95% confidence level, it leads on . However, a DTS required a large computing time and the realization of simulations is often not feasible. This paper follows an approximate approach, proposed by the same author: a relatively small value of simulations is performed. The average and the standard deviation of the outputs are computed to characterize the distribution of the model outputs.

The most commonly used outputs in UP, are some global results as the building consumption during one year [START_REF] Tian | A review of uncertainty analysis in building energy assessment[END_REF], the maximal heating load [START_REF] Tahmasebi | A systematic assement of the sensitivity of building performance simulation results with regard to occupancy related input assumptions[END_REF] or the cumulative distribution of indoor temperature [START_REF] Mahnameh | Consideration of Inhabitants Diversity in Building Performance Simulation: Does It Matter?[END_REF]). Yet, the temperatures of the thermal zones are also interesting outputs of DTS. Recent works such as [START_REF] Maykot | A field study about gender and thermal comfort temperatures in office buildings[END_REF] or [START_REF] Aghniaey | Thermal comfort evaluation in campus classrooms during room temperature adjustment corresponding to demand response[END_REF] used the measurement of the inside air temperature for a thermal comfort analysis of the building. In an experimental process, the temperature can be a key factor to identify parameters related to the building energy performance [START_REF] Roels | On site characterisation of the overall heat loss coefficient: Comparison of different assessment methods by a blind validation exercise on a round robin test box[END_REF][START_REF] Ajib | Data-driven building thermal modeling using system identification for hybrid systems[END_REF]. That is why this contribution focuses on these outputs only. At the end of the UP, several methods exist to analyze the model outputs variability. The standard deviation and the average of the observations of each output can be computed [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF]. However, if the number of output is too large, it will be difficult to summarize these results [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF]. For this purpose, [START_REF] Campbell | Sensitivity analysis when model outputs are functions[END_REF] have reduced the dimension of the outputs space with a principal component analysis. When just a few axis can be retained, this technique allows to interpret easily the UP results. When more axis has to be kept, the interpretation is much more difficult.

Method Overview

In the operation phase, a building monitoring rarely gives a very precise knowledge about the internal heat gains. The lack of information is supposed to be related to:  The measurement of the total internal heat gains of the building ;  The spatial repartion of the internal heat gains between the thermal zones of the DTS model ;  The time dependence of spatial repartition. In order to evaluate the effect of these different uncertainty sources, three UPs are applied to a DTS building model. They take into account independently the uncertainty sources previously mentioned. At the end of one UP, a large number of temperature time series are simulated. A pre-selection step allows to focus the analysis on a part of the outputs whose variances are representative of all the outputs. This step operates in two times:  First, by finding temporal patterns in the variance of the average temperature of the building; it allows to select shorter periods representative of the entire period of simulation;  Secondly, the influence of each uncertainty source on the temperature of each thermal zone is assessed by calculating cumulative variances. From these data, the thermal zones are grouped with a clustering algorithm and one thermal zone is selected as a representative of each group. Finally, the average and the standard deviation of the distribution of the pre-selected outputs are analyzed.

Simulation model

The model is a multi-zone model. It divides the building into thermal zones. Each sub-volume is depicted as a uniform temperature and air pressure volume. Each room of the building is considered as a zone of the model. The simulations were performed with (TRNSYS 17, 2010) software and CONTAM software [START_REF] Dols | CONTAM User Guide and Program Documentation Version 3[END_REF] to couple heat and air transfer. These tools were chosen because they allow the prediction of the dynamic of the air temperature of several thermal zones. During the heated period, a model of an ideal heating system controls each room temperature. This model perfectly adapts the indoor temperature to the temperature set point thanks to a direct acting-control. The building performance monitoring (as described in 'Case study' part) provides information for the parameterization and the definition of its inputs. The internal heat gains of the model are stored in an internal heat gains matrix defined with:

[ ( ) ( ) ( ) ( ) ( ) ( ) ]
(1) With , the simulation timestep and , the period of simulation. The simulation discretizes a one year period into = 8760 hours. is the index of the modelled thermal zone and is the number of thermal zones. The internal heat gain of the entire building, quoted ( ) is the sum of the columns of . The part of the internal heat gains assigned to the thermal zone is called an internal heat gain ratio (noticed ( )). In the definition of , assumptions have to be made about the time variation of these ratios (fixed ratios or time dependent ratios). The simulated temperatures are stored in matrix where each element ( ) corresponds to the temperature of the thermal zone at the time .

[ (

(2)

Uncertainty propagation

First, the probability distribution of the internal heat gains parameters has to be quantified according to the uncertainty source. [START_REF] Rivalin | A comparison of methods for uncertainty and sensitivity analysis applied to the energy performance of new commercial buildings[END_REF] advice to choose normal distribution for measured parameters. For a non-measured parameter, if an expert can suggest a minimum and a maximum value, it will follow a uniform distribution.

The 1 st UP is related to the measurement uncertainty of ( ).

( ) is estimated with the measurement of the electricity consumption of the entire building and an occupancy schedule built with the building manager. Then, the th observation of ( ) affected by the measurement uncertainty, is written ( ), and follows:

( ) ( ) ( ) (3) with (
) where is the mean and is the standard deviation of a normal distribution. These parameters have been quantified with a measurement uncertainty analysis by considering the uncertainty on the part of electricity consumption reused as a building heat flow, the lack of knowledge about the number of people and their metabolism [START_REF] Titikpina | Optimisation des plans de mesure dans un contexte de garantie de performance énergétique des bâtiments[END_REF]. The generation of random draws of is stored in the 1 st input sample called . The 2 nd UP is related to the uncertainty about the spatial repartition under the assumption that internal heat gain ratios are fixed in time. These ratios are defined with expert judgment according to the functionality and the available information about the number of occupants in each thermal zone. To take into account the uncertainty on that method, each ratio is supposed to follow ( ) with and the parameters of a uniform distribution. The generation of random draws of (the vector of all the fixed ratios) is the 2 nd input sample called . Finally, the 3 rd UP is related to the uncertainty about the spatial repartition under the assumption that internal heat gain ratios are time dependent. However, no known method allows modelling the time evolution of the internal heat gain ratios. It can only be pointed out that some authors proposed stochastic approaches or random walk approaches for modelling the random nature of the occupant presence [START_REF] Ahn | Different Occupant Modeling Approaches for Building Energy Prediction[END_REF]. In a first approach, a purely random process is used for modelling ( ). Each ratio at time follows ( ) ( ). The uniform distribution parameters are the same as the case of . The generation of random draws of ( ) (the matrix of the random process of all the ratios) is stored in the 3 rd input sample .

The output of the simulation of the UP is noticed . Each element ( ) of is the temperature of the thermal zone at time for the simulation . Then, the estimated variance of the response ( ) is noted ( ) and is computed with:

( ) ∑ ( ( ) ( )) (4)
where ( ) is the estimated average of ( ) which is equal to:

( ) ∑ ( ( )) (5)
( ) quantifies the variability of each temperature caused by the uncertainty source. Outputs pre-selection First, the average temperature of the building ( ) is computed with:

( ) ∑ ( ( )) (6)
By combining equations ( 6) and ( 4), the variance of the average temperature ( ) is obtained. The variations over the time of ( ) show the time evolution of the influence of the uncertainties sources. Then, some temporal patterns of ( ) may exist and are representative of the whole period. Then, the next step will focus on these shorter periods. The second step is the selection of representative thermal zones. First, the influence of each uncertainty source on each zone is assessed by the calculation of cumulative variance [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] on pre-selected periods :

∑ ( ) (7)
where (resp. ) is the beginning (resp. the end) of a shorter period.

indicates the influence of one uncertainty source on the temperature of the zone during the period between the interval . Then, an exponent is added to in order to specify the period and the source of uncertainty. Thus, each thermal zone is characterized by the cumulative variances computed under the different conditions previously described. The cumulative variances of each thermal zone are stored in the cumulative variances matrix :

[ ] (8)
Each zone is associated with the row of . Each column of is normalized. The Euclidean distance matrix between each row is computed. Then, a complete linkage method clusters the thermal zones according to a hierarchical scheme. This data analysis was computed with (R Development Core Team, 2018). Finally, for each cluster, one thermal zone is chosen as a representative of the others elements of the same group.

Uncertainty analysis

The distribution of each pre-selected output is analyzed by representing their confidence interval ( ) with a 95% confidence level. Under the assumption that ( ) is normally distributed follows:

[ ( ) √ ( ) ( ) √ ( ) ] (9)

Case study

The case study is an office building (Figure 1). It was constructed from 2010 to 2012 near to Angers (France).

It has two floors. Two parts divide the building according to the usage. One part is the office (853m²) and the second part is a workshop (703m²). The building complies with French Thermal Regulation RT2005 that was in effect during its design. However, the performance of the 'office' part is very close to the criteria of the Low Consumption Building -Effinergie label (BBC -Effinergie, 2007). The 'office' part was modelled with TRNSYS software [START_REF] Titikpina | Optimisation des plans de mesure dans un contexte de garantie de performance énergétique des bâtiments[END_REF]. The model divides the building into thermal zones that correspond to the same room division of the 'office' part of the building. weather is monitored. These data are used as inputs of the DTS. Moreover, electricity meters record each hour the building electricity consumption. Combined with the knowledge of an hourly occupancy schedule of the building, these data are used to define the internal heat gains of the entire building ( ( )). The ratios are defined relatively to each another, according to:  The functionality of the thermal zone: for instance, an office room is fitted with personal computer and office worker can be present in this zone. The ratio will be more important in this room than another room like a corridor;  The specific knowledge about the real occupancy of the room: for instance, during the monitored period, BUR006 was not occupied. This information was taken into account to affect a lowest ratio to this room. [START_REF] Titikpina | Optimisation des plans de mesure dans un contexte de garantie de performance énergétique des bâtiments[END_REF] applies an uncertainty measurement analysis to this case study that shows ( ). The determination of the parameters of the uniform distribution of is made with an expert judgment seeking to conserve the initial information about the spatial repartition.

Results & Discussion

Outputs pre-selection Figure 2 shows the variance of the average temperature of the building ( ( )). First, it can be noticed that the variance estimated with is larger than and . The uncertainty on the entire internal heat gains leads to a larger variability than the uncertainty related to the spatial distribution (with fixed or time varying ratios).

( ) follows week variation. It may be caused by the week variation of the entire internal heat gains. On the other hand, it can be noticed that the model response to is strongly influenced by the non-heated period. These observations lead to define two shorter periods: « summer» and « winter » as shown in Figure 2. 2016)) of the cumulative variances matrix ( ). Each row refers to a thermal zone . Each column refers to a particular uncertainty propagation (defined with , or ) and a specific period ('SUM' stands for the summer period and 'WIN' stands for the winter period). The clustering results are represented on the Figure 3. The rows of were reorganized in order to position the zones close to the ones with similar cumulative variances. conditions. Figure 4 shows that the selected thermal zone BUR_101 seems close to the selected thermal zone of G1. There are few distinctive elements between these two objects;  G5 regroups the last 10 zones. Globally, they are the less impacted zones for all the conditions. Now, a sufficient number of representative outputs are selected. The next subsection will describe the uncertainty analysis of the chosen ones.

Figure 4 : Normalized cumulative variances according to the selected thermal zones Uncertainty analysis results G1 gathers a majority of thermal zones relied on individual room office (actually, the 'BUR' prefix refers to a room office). Indeed, that kind of rooms is often used at different levels of evaluation of the building thermal performance, and thus, this is a good subject of interest. BUR_106 (G1) is an individual office (12m²) on the 1 st floor. It is located between two others individual offices and is accessible from the main corridor of the respective floor. This room has a north facing facade. One person occupied this office during the data-monitoring period. Figure 5, Figure 6 and Figure 7 show the distribution of the temperature of BUR_106 simulated with different input uncertainty sources and during different periods. The black line is the mean value. The grey area covers the confidence interval computed with (9). Figure 5 shows that the distribution width is affected by the uncertainty of the entire heat gains when the room temperature is not controlled by the heating system. However, during the summer period, the distribution of the temperature around the average value is fixed throughout the all period (Figure 6). When the influence of the spatial distribution is tested, the distributions are tightening up (Figure 7). The uncertainty related to the distribution has less impact than the uncertainty of the entire heat gains. The thermal zone SANIT_HOM_R1 belongs to (G2).

The distributions are different to BUR_106, especially for the winter period. As this is a non-heated room, the temperature is influenced at each time step of the winter period (Figure 9). SAL_ESSAI belongs to the group (G3). This zone is identified as the most sensitive zone to the variation of the spatial distribution during the summer period. Figure 8 shows the temperature distribution of this zone obtained with M2. The temperature seems mainly influenced by the spatial repartition during the workdays. G4 gathers the remaining building offices. The distributions of BUR_101 (G4) are close to the outputs distributions of BUR_106 although the clustering methods lead to two different groups.

between two non-heated zones

The thermal zone SANIT_FEM_R1 belongs to (G5). As SANIT_HOM_R1, it is a non-heated room too, but the model response is less impacted during the winter period (Figure 9).

Average of standard deviation

The simulated temperatures could be useful in order to analyze the thermal comfort of every zone of that building. However, the UPs show the simulated temperatures could be sensitive to the accuracy of the internal heat gains. The influence also depends on: the time period, the thermal zone and the uncertainty source. Figure 10 represents the average (in each period) of the standard deviation deducted from the cumulative variances ( 7). The highest value of the average of standard deviation is related to the uncertainty on the entire heat gains during the summer period. It is observed for BUR_106 and is equal to 0.4°C. During the winter period, the hisghest influence is obtained for a non-heated zone (SANIT_HOM_R1) and is equal to 0.2°C. The accuracy of a thermal comfort analysis based on that DTS could be improved by enhancing the accuracy of the measurement of the entire heat gains. The uncertainty on spatial distribution could lead on an average of standard deviation equal to 0.2°C (SAL_ESSAI during the summer period). However, for the same thermal room and for the same period, in considering that the internal heat gains ratios change randomly in time, the average of standard deviation reaches a value of 0.07°C.

Conclusion

This paper has shown a method to investigate the effect on the building thermal behaviour of different internal heat gains uncertainties (the heat gains of the entire building and their spatial distribution). Input samples built with the Monte-Carlo approach have been used to estimate the variance of the modelled temperatures of 39 thermal zones over a period of one year. This paper has proposed a two steps method to select the most interesting part of theses outputs.

An office building has been used as a case study. In that case, the uncertainty source that has the most influence was generally the uncertainty on the measurement of th SANIT_HOM_R1 SANIT_FEM_R1 entire heat gains. However, the spatial repartition has had a larger influence on the thermal zone temperature with some specific conditions, which are in the presented case: a winter period, a non-heated zone and a spatial repartition fixed in time.

In the case study, the measurement uncertainty of the entire internal heat gains was the most influential parameter. It was mainly influent during the summer period of the simulation. During this period, the average of standard deviation could reach 0.4°C. It was shown that, if the internal heat gains ratio follows a random process, then, the uncertainty related to the spatial distribution has less impact. This implication can reduce the average of standard deviation from 0.2°C to 0.07°C

The variance of the model outputs depends on the model parameterization and the implementation of the uncertainty propagation. For instance, the cumulative variances of the response are different for a non-heated or a heated zone during the winter period. The identification of more complex features could be a satisfying improvement of this methodology.

To complete this work, the interaction between the entire heat gains uncertainty and the spatial distribution uncertainty should be reviewed. On the other hand, this work could be improved by observing relationships between the input sample factors and the model responses with a sensitivity analysis.

Figure 1 :

 1 Figure 1: Picture of the case study office building The simulation computes the building thermal behavior from the 2015/03/01 to the 2016/02/28. At the end of each simulation, the 39 temperature time series are recorded. During the simulation, a heated period is allocated between 01/03/2015-17/05/2015 and 27/09/2015-28/02/2016. The real not-heated rooms of the building are configured in the model. The local

Figure 2 :

 2 Figure 2: Time evolution of ( ) according to the uncertainty propagation ( ) just gives a global view of the model responses at the building scale. The analysis of the cumulative variances at the thermal zone scale is used to choose representative thermal zone. Figure 3 is a heatmap representation (realized with d3heatmap, a R package

Figure 3 :

 3 Figure 3: Heatmap and cluster ofgreen cells: high level of cumulative variances, purple cells: low level of cumulative variances It seems relevant to divide the 39 zones into 5 different groups. A representative of each group is arbitrarily selected: BUR_106 for G1, SANIT_HOM_R1 for G2, SAL_ESSAI for G3, BUR_101 for G4, SANIT_FEM_R1 for G5. Figure 4 shows the cumulative variances of each of them. It helps to understand the properties of each group:  G1 regroups 19 zones. It regroups the most influenced zones by the variation of and during the summer period. They are moderately impacted by the others conditions;  G2 regroups 1 zone : SANIT_HOM_R1. This zone is very affected by the variation of and during the winter period;  G3 regroups 2 zones. They are the most influenced zones by the variation of and during both periods. They are moderately affected by ;  G4 regroups 7 zones. They are affected by the variation of during the summer period. Moreover, they are the less influenced zones by the other

Figure 5 :

 5 Figure 5: Responses of BUR106 to M1 -Winter period

Figure 8 :

 8 Figure 8 : Responses of SAL_ESSAI to M2 -summer period

Figure 10 :

 10 Figure 10 : Average of standard deviation according to the uncertainty sources and the simulated period