Jessica Rousseau
email: jessica.rousseau@hmg.inpg.fr

Emmanuel Frangin

Philippe Marin

Laurent Daudeville

SOME ASPECTS OF DE/FE COMBINED METHOD: MODEL IDENTIFICATION AND COMBINED METHOD

Keywords: Discrete Elements, Finite Elements, Concrete, Impacts

INTRODUCTION

This study deals with impacts on reinforced concrete structures. Its aim is to predict local damage in the vicinity of the impact as well as the global response of the structure. The Discrete Element Method (DEM) seems particularly convenient for the modeling of fractures. The chosen Discrete Elements (DE) are rigid spheres of different sizes. The identification of the DEM material parameters from macroscopic data (Young's modulus, compressive and tensile strength, fracture energy…) will be presented in the first part.

The modeling of a large structure by means of DE may lead to prohibitive calculation durations. A refinement is required in the vicinity of the impact whereas the structure can be modelled by means of a coarse FE mesh far away from the impact area, where the material behaves elastically. Thus a coupled discrete-finite element approach is proposed: the impact zone is modeled by means of DE and elastic FE are used for the rest of the structure. The proposed approach will be described. Finally, an application of this coupling method will be shown by means of the response of a concrete slab impacted by a rocky block.

IDENTIFICATION OF THE DISCRETE ELEMENT MODEL PARAMETERS

Discrete Element Model

The used Discrete Element method is fully described in Hentz [START_REF] Hentz | Identification and validation of a discrete element model for concrete[END_REF] . Interactions between two spheres are defined with normal K n and tangential K s stiffness characterizing the elastic behavior of concrete. "Micro-macro" relations (1) are supposed in order to compute the stiffness Kn and Ks between elements a and b from Young's Modulus and Poisson's ratio. Parameters α,β,γ of these relations were identified for a given average number of interactions per discrete element.

A modified Mohr-Coulomb criterion associated with softening is used to model the nonlinear behavior of the material. Local parameters, as local tensile strength T, cohesion Co and softening factor ζ need to be identified from global parameters as compressive and tensile strength σ c and σ t and fracture energy G f . A classical friction constitutive behavior is used between elements in contact.

, 1 K s a b K i n i t n n K s i n t K n K s K n K s K n D E K S β γ α ν α  +  =  +    -  =  +  or ( ) , 1 1 (1 ) 1 1 i n t n a b i n i t s n S K E D K K α β ν γ α ν α ν ν +  =  + + -   -  =   + (1)

Parameters identification of linear elastic constitutive behavior

Simulations of a compressive test with a DE sample are realized with different values of K s /K n . Numerical values of E and ν are computed and values of α, β,γ are then identified with (1) in order to fit numerical points of E and ν versus K s /K n .

Tests have been carried out with 15 different samples (10 of them have identical sizes but different disorder and the 5 others have different geometries). The main assumption of this approach is that samples are representative of modelled concrete material in further simulations of impact on concrete structures. Samples have to be isotropic. With the identified values of α, β,γ parameters, average uncertainty is about 2% (maximal error is 7%) for Young's modulus and about 3% (maximal error is 6%) for the Poisson's ratio. In this part, the point is to identify local parameters to model macroscopic values as compressive strength σ c , tensile strength σ t and fracture energy G f . Only local tensile strength T, cohesion Co and the softening factor ζ will vary in this study, friction parameters will be fixed values.

The aim is to realize traction and compression simulations in order to simulate tests performed in 3S-R laboratory (Gabet 2 ). We identified parameters according to Hentz's procedure. In Figure 1, we can observe the result of the traction simulation and the compressive tests on Figure 2.

IMPACT ON STRUCTURE WITH COMBINED DE/FE METHOD

Description of the method

This third part concerns the coupling between continuum and discrete domains. Discrete element method is convenient for modeling discontinuities and fracture but it may lead to prohibitive durations of calculation with a large structure. During an impact, though the impact zone needs a refined mesh, the main part of the structure far away from the impact zone can be modeled by means of coarse finite elements. We propose to couple discrete elements in the vicinity of the impact with elastic finite elements in the rest of the structure. It is chosen to use a bridging domain in order to couple DE and FE degrees of freedom.

Many methods have been already developed for multi-scaled problems. Belytchko 3 proposed to use a bridging domain where the Hamiltonian is a linear combination of discrete and finite ones. Equation (3) shows DE and FE Hamiltonian are weighted by a parameter α equals to 1 in the FE domain and equals to 0 in the DE domain.

( 1)

FE DE H H H α α = + - (3)
We link in the bridging domain DE displacements d r with FE displacements u r . DE rotations ω r are linked with the antisymmetric part of the FE displacement gradient . The solution is obtained by minimizing the Hamiltonien with introduction of the kinematics constrains by means of two Lagrange multipliers linked to displacement and rotation constrains. Details of the method are explained in Frangin [START_REF] Frangin | On the use of combined finite/discrete element method for impacted concrete structures[END_REF] .

Due to the variation of discretization size between both models, we use a process based on the relaxation of Lagrange multipliers or on damping parameter to reduce the high wave reflection (Frangin 4 ).

Impact on a concrete structure

To quantify the efficiency of this method, we propose to compare damage and displacement response of a concrete slab under a rocky impact (Figure 3) between combined FE/DE method and the full DE method (Figure 4). Displacements and damage of the slab are very close at each part of the slab. The combined method is 10 times faster then the fully DE method without a multi-time step algorithm. Some computation on a reinforced concrete structure will be presented. method and coupled method (blue)

Influence of the reinforcement

Main modification induced by the modeling of reinforcement concerns the random process where DE cannot move easily around reinforcement and the porosity of the media is higher near the reinforcement than in the concrete. Two different interaction radii are used to keep a number of contact equals to 12 for concrete DE and for steel DE. The FE domain and the bridging domain keep an elastic behaviour and so we do not take into account the reinforcement in the FE domain. The bridging domain is a way to regularize the displacement and there is no localization due to the reinforcement in this area.

Figure 1 :

 1 Figure 1 : Uniaxial Traction : Stress-Strain curve Figure 2 : Uniaxial Compression : Stress-Strain curve

Figure 3 :

 3 Figure 3: The discrete model Figure 4: Comparison of displacement between DE (red)

CONCLUSIONS

-

The identification of parameters for concrete behavior gives very accurate results. First, the elastic behavior seems to be reliable. And secondly, uniaxial tensile and compressive behaviors can be easily computed. -This identification process combined with an efficient coupled DE/FE method give solution to simulate large RC structure under impact.


Xingheng Liu

Yann Dijoux
email: yann.dijoux@utt.fr

Jørn Vatn

H Åkon Toftaker

PERFORMANCE OF PROGNOSIS INDICATORS FOR SUPERIMPOSED RENEWAL PROCESSES

Keywords: applied probability, reliability theory, simulation, stochastic modelling

The paper deals with prognosis estimation for industrial systems in a series configuration, modeled by superimposed renewal processes (SRP), when the cause of failures is not available. In the presence of missing information, an SRP is commonly approximated by a Poisson process or a virtual age model. The performance of the approximations was assessed in the ideal configuration where all parameters of the models are known. The current article adopts a practitioner's perspective by assuming that the parameters of the models are unknown and must be estimated. In addition to inference procedures, the assessment of the prognosis indicators, such as the remaining useful life, is discussed. Finally, we investigate a fleet of infrastructure components of the Norwegian railway n e t w o r ko p e r a t e db yB a n eN O R .

INTRODUCTION

In this paper, the superposition of independent renewal processes (SRP) is discussed in the context of reliability engineering [START_REF] Song | An integrated method for estimation with superimposed failure data[END_REF][START_REF] Wu | A failure process model with the exponential smoothing of intensity functions[END_REF][START_REF] Zhang | Estimating a parametric component lifetime distribution from a collection of superimposed renewal processes[END_REF] where it can be referred to as partial repair model [START_REF] Bedford | Competing risks in reliability[END_REF][START_REF] Bedford | The identifiability problem for repairable systems subject to competing risks[END_REF]. A maintainable system under partial repair is a series system such that after a failure, the failed component is automatically identified and instantaneously replaced by a new one of the same technology (as good as new). The non-failed components are assumed to be minimally maintained (as bad as old) which implies that their respective failure rate is the same before and after the maintenance activity. At a system level, this partial repair is also characterized as imperfect repair [START_REF] De Toledo | ARA and ARI imperfect repair models: estimation, goodness-of-fit and reliability prediction[END_REF][START_REF] Tanwar | Imperfect repair modeling using Kijima type generalized renewal process[END_REF]. The components are assumed to be statistically independent and can represent either physical units or virtual ones [START_REF] Nafisah | Virtual series-system models of imperfect repair[END_REF] associated with an independent failure mode.

Maintenance records are often incomplete, and the information on the identity of the failed components is not necessarily available. Consequently, the observations are commonly reduced to a pooled output [START_REF] Cox | On the superposition of renewal processes[END_REF], consisting of the failure times at a system level. Assessing the health of the system such as its overall aging and the maintenance efficiency usually starts off by estimating model parameters. When the number of components is known and all components are identical, inference procedures have been carried out by Zhang et al. [START_REF] Zhang | Estimating a parametric component lifetime distribution from a collection of superimposed renewal processes[END_REF] directly from an SRP using computational partitioning. Nevertheless, when the number of pooled events is limited, or when the components are different from each other, direct inference methods are hard to implement, which inspires diverse approximation approaches.

Alternately, an SRP is approximated by a homogeneous Poisson process [START_REF] Proschan | Mathematical theory of reliability[END_REF], a renewal process [START_REF] Whitt | Approximating a point process by a renewal process, I: two basic methods[END_REF], an imperfect maintenance model [START_REF] Liu | On approximation of superposition of renewal process[END_REF][START_REF] Wu | A failure process model with the exponential smoothing of intensity functions[END_REF] or a copula model [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF]. The benefits of preferring imperfect maintenance models or copulas to a renewal process, when the parameters of the approximated SRP are known, have been highlighted in [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF]: these models account for the correlation between consecutive inter-occurrence times. In this paper, we shall evaluate the above-mentioned approximation approaches when the model parameters are estimated from the pooled output.

The contribution of the article is multi-fold. First, an innovative procedure to evaluate the performance of five different approximation approaches is presented based on density level sets. Second, the estimation procedures implemented for each approximating model of an SRP are derived and compared. The third contribution consists of taking into account periodic preventive maintenance within the pooled outputs. Condition-based maintenance requires more complex modeling and is beyond the scope of the paper. Finally, we emphasize the selection of particular approximation methods along with specific inference procedures on a real data set from Norwegian railway signaling systems.

The rest of the paper is structured as follows. In Section 2, the SRP notations are introduced before presenting a new performance-measuring procedure based on level-set. Section 3 addresses the inference procedures and evaluates the performance of the approximations. An application to real data is proposed in Section 4. Most systems in the Norwegian railway signaling network operated by Bane NOR can be assumed to be series systems under partial repair with both corrective and potential preventive maintenances. The benefits of the aforementioned approximations are highlighted through the modeling and estimation procedures based on the pooled output of signaling systems and relevant prognosis indicators.

APPROXIMATION OF A SRP

Notation

Consider a maintainable series system under partial repair. A corrective maintenance (CM) is carried out immediately after a failure and consists of replacing the failed component while maintaining the others minimally. Periodic preventive maintenance (PM) is performed at ∆, 2∆, 3∆ ...,a n d∆=∞ if no PM is implemented. It is also assumed that at t =0, the system is in its stationary regime, that is, has been functioning for a long while. The identity of the components causing the failures is assumed to be unknown. Observations can be modeled by a colored point process: failure times are denoted by {T i } i≥1 , with inter-arrival times {X i } i≥1 . The indicator of the maintenance types are {δ i } i≥1 with δ i =0 for CM and δ i =1forPM.Figure 1 presents an hypothetical trajectory of a maintenance process for a dwarf signal, wherein the lamps are in series configuration: the four lights function independently and are subject only to CM, and no identification regarding the renewed bulbs is available. It is, therefore, reasonable to model the pooled failure times {T i } i≥1 by an SRP.

Reminder of the Current SRP Approximations

SRP are usually approximated by a standard renewal process [START_REF] Torab | On approximate renewal models for the superposition of renewal processes[END_REF][START_REF] Whitt | Approximating a point process by a renewal process, I: two basic methods[END_REF]. In a first study [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF], we have highlighted the benefits of using virtual ages or copulas to approximate an SRP:

• A virtual age model [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] is an imperfect maintenance model which assumes that after maintenance at time T i , the system behaves as a new and unmaintained one of age A i . The virtual age model that approximates plausibly an SRP when components are homogeneous is the arithmetic reduction of age with infinite memory, ARA ∞ [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF][START_REF] Nguyen | Analytical properties of an imperfect repair model and application in preventive maintenance scheduling[END_REF], which assumes that the virtual age of the system after a maintenance is reduced proportionally:

A i =(1-ρ)(A i-1 + X i ) ,ρ ∈ [0, 1]. ( 1
)
ARA ∞ is a stable process [START_REF] Liu | Steady-state imperfect repair models[END_REF] in the sense that its intervals converge in distribution. Like in an SRP, the failures in an ARA ∞ model arrive at a constant rate, asymptotically.

• A copula model [START_REF] Joe | Dependence modeling with copulas[END_REF][START_REF] Sklar | Random variables, distribution functions, and copulas: a personal look backward and forward[END_REF] is a convenient approach to characterize a joint distribution by defining the marginal distributions and a copula function. The copula associated with the observation of two consecutive inter-arrival times of an SRP in its steady regime (X i ,X i+1 ) has been employed to approximate an SRP. This modeling is automatically more efficient than a renewal process as it characterize the distribution of X i+1 given X i . This model is denoted by inter-arrival-time copula of memory 1 (IAT 1 ).

When the parameters of the approximated SRP are known, ARA ∞ and IAT 1 outperform the renewal process approximation for giving more accurate estimation of the remaining useful life (RUL) [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF], with the ARA ∞ being marginally better than IAT 1 .I n addition, both ARA ∞ and IAT 1 capture the negative dependency between two consecutive inter-arrival times when all the components have an increasing failure rate (IFR) [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF]. The correspondence between SRP and ARA ∞ can be emphasized further using one of the strongest measures of negative dependence called reverse regular of order 2 (RR 2 ) [START_REF] Balakrishnan | Continuous bivariate distributions[END_REF]. The proof is presented in Appendix A.

Proposition 2. 1: Two successive intervals in an SRP, formed by components with increasing failure rates, are RR 2 -dependent. Similarly, two successive inter-arrival times of an ARA ∞ process with increasing baseline failure intensity are RR 2 -dependent.

A Level-Set

=[T i + a ′ ,T i + a ′ + b -a]o f amplitude b -a
can be determined which maximizes the confidence set of the approximating model. The probability that T i+1 belongs to I is naturally 1 -α, and the probability that T i+1 belongs to J is automatically less than 1 -α. The quality of the approximating model lies in how the latter probability is close to 1 -α. This probability p α can be estimated empirically from a unique and sufficiently long trajectory by assessing the proportion of times when the prediction has been correct. Furthermore, the overall performance of the approximated model can be synthesized by its Gini index G =2 1 0 (1 -α -p α )dα.IfG =0, the approximating model is indistinguishable from the actual model, and G = 1 is the worst-case scenario where the approximating model is totally inconsistent with SRP.

Five approximated models are considered, and their parameters are assumed to be known in this section only:

HPP. The homogeneous Poisson process (HPP) with a rate is equal to the asymptotic rate of the SRP. This model is the simplest SRP approximation and is valid when the number of sources in the SRP tends to infinity [START_REF] Drenick | The failure law of complex equipment[END_REF].

SIM. The stationary interval method (SIM; [START_REF] Torab | On approximate renewal models for the superposition of renewal processes[END_REF]) corresponds to a renewal process with a generic distribution equal to the asymptotic distribution of the SRP.

IAT 1 . Inter-arrival-time copula of memory 1, as reminded in Section 2.2. ARA ∞ . ARA ∞ process, as reminded in Section 2.2.

BP. The Brown-Proschan model (BP; [START_REF] Brown | Imperfect repair[END_REF]). Belonging to imperfect maintenance models, BP assumes that maintenance is either perfect with probability p or imperfect with probability 1 -p. Like an SRP, BP has a stationary regime [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF]. It is conceivable that SRP with heterogeneous components are relatively better approximated by BP, wherein the maintenance is heterogeneous. The parameters chosen for the BP model are the shape and scale parameters of the initial Weibull failure rate and the parameter p which has the best goodness-of-fit (GOF) with the initial SRP.

Five SRP configurations are studied as follows: I. A system formed of three homogeneous components with a low wear-out rate. The failure rate of each component follows a Weibull distribution W(η, β) with a shape parameter β =1.5 and a scale parameter η =1. V. A system formed of six relatively homogeneous components with a moderate wearout rate. The components have a Weibull lifetime distributions: three are W(η = 1,β =2.5) and three are W(η =2,β =2.5).

For each model and each configuration, the level set p α is computed and compared to the initial confidence level 100(1 -α)% through the identity line. Figures 2 and3 present the results with a fast wear-out for the configurations II and IV, respectively. The figures for other configurations are not presented as the plots are almost indistinguishable from the identity line. In addition, each Gini coefficient is provided in Table 1.

The analysis of the level sets and Gini index completes and refines the initial analysis from [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF]:

• For each model, the Gini index is relatively low, between 0.002 and 0.2, which indicates that each model provides an efficient approximation of an SRP. The Gini index are particularly low when the aging of the components are moderate or when the number of components is important.

The model under ARA ∞ assumption almost consistently provides the best index, regardless of the heterogeneity of the components and the aging rates of the components. • The IAT 1 and BP models present satisfactory but more contrasting results. The BP approximation is comparatively much less efficient when the components are aging fast (configurations II and IV). The IAT 1 approximation outperforms the other models in the heterogeneous case with fast aging rates (configuration IV) but provides, otherwise, similar results to the HPP approximations.

• The two static models HPP and SIM presents decent and very similar results. When the aging rate is moderate, their efficiency has the same order of magnitude as the other adaptive models, but they are clearly outperformed by ARA ∞ when the wear-out rate is fast.

• All the models are extremely efficient when the number of components is large, particularly the HPP which is the theoretical limiting case when the number of components tends to infinity [START_REF] Drenick | The failure law of complex equipment[END_REF].

PERFORMANCE OF THE APPROXIMATIONS WITH UNKNOWN PARAME-TERS

The parameters of the models are commonly not known and need to be estimated first from a relatively small data set. The performance of the inference procedures might put the previous analysis into perspective.

Configurations

The performances of the SRP approximations are discussed for the five models: HPP, SIM, IAT 1 ,ARA ∞ and BP. The number of observations is assumed to be larger than 10 without exceeding 200, and the system is assumed to be in its stationary regime at the beginning of the observation. Three configurations previously studied in Section 2.3 are considered, two with CM only and one with CM and periodic PM:

(A) The configuration II with three homogeneous components and a relatively fast wearout (W(η =1,β =3.5)) and CM only.

(B) The configuration IV with three heterogeneous components and a relatively fast wear-out (W(η i = {1, 2, 10},β =3.5)) and CM only. (C) The configuration II with three homogeneous components and a relatively fast wearout with both CM and PM. The periodicity of the preventive maintenance policy is ∆ = 0.425. The periodicity has been chosen so that the proportion of PM in the actual SRP model is 70%. The maintenance efficiency is assumed to be imperfect but of great quality on each component. Therefore, it has been opted to assume that PM efficiencies are ARA ∞ with an improvement factor ρ =0.7. The virtual age of each component after a preventive maintenance action is 30% of its virtual age just before the PM.

These three configurations are not exhaustive but are quite representative of the behavior of an SRP with or without homogeneity. The PM policy in (C), which corresponds to the case where imperfect PM are carried out simultaneously on all the components, is investigated because it is easy to implement and widely applied in practice.

Measures of Performance

The objective is to quantify the efficiency of the approximation approaches when the number of observations is limited. Two predictors, the mean lifetime E[X ∞ ] and the Pearson's correlation coefficient Corr between two successive inter-failure times, have been considered. They are assumed to be known for the SRP models, either computed theoretically considering CM only or obtained based on Monte Carlo simulations at any given precision when both CM and periodic PM are implemented. Given a sample size N of an SRP configuration, the mean squared error (MSE) of the two predictors E[X ∞ ] and Corr are derived empirically for the five approximated models based on 5,000 replicates of histories. A general structure of the simulation procedure is presented in Algorithm 1.

Inference Procedures

Additional information is provided on the estimation methods for the five models corresponding to lines 7 and 8 in Algorithm 1. In the case of no PM, the observations consist of the inter-failure times {X i } i=1...N or equivalently the maintenance times {T i } i=1...N .Inthe case of both PM and CM, the observations additionally comprise the maintenance types {δ i } i=1...N .

Homogeneous Poisson process

The failure rate λ in an HPP can be assessed by MLE: for i=1:5 //ith model (HPP, SIM, IAT 1 ,A R A ∞ ,B P )

λ = n - n i=1 δ i T n . (2)
7:
Assess the parameters of the ith model from Hist

8:

Estimate the predictors by e i (k) and r i (k)

9:
end for 10:

end for 11:

for i=1:5

12:

Compute the empirical MSE of each measure for a trajectory of length N and the approximated model i:

13: MSE E (i, N )=1/5000 5000 k=1 ( e i (k) -e) 2 14: MSE Corr (i, N )=1/5000 5000 k=1 ( r i (k) -r) 2 15:
end for 16: end for The first predictor, E[X ∞ ], equals 1/ λ. The second predictor, Corr, is 0 since the intervals are independent in an HPP.

Stationary interval method

The SIM method approximates an SRP by a renewal process whose interval is distributed identically as X ∞ . This distribution can be roughly estimated by the empirical survival function in case of CM only or by the Kaplan-Meier estimator function in the presence of CM and PM.

To evaluate the accuracy of SIM approximation, E[X ∞ ] is computed by numerically integrating a piecewise-linear version of the estimated survival function, and the correlation is 0 given that successive inter-arrival times are independent for a renewal process. Within the course of the inference procedures, additional non-parametric estimations of the reliability function have been tested using different kernels and several bandwidths, but the piecewise-linear version of the Kaplan-Meier has been opted as it consistently provided the best results.

In the presence of PM, the Kaplan-Meier estimator is not defined beyond the largest observation if the longest interval is censored, which makes it impossible to compute the expected value. We adopt therefore the common approach, that is change the largest observation to a death time if it is censored [START_REF] Efron | Assessing the accuracy of the maximum likelihood estimator: observed versus expected fisher information[END_REF].

Inter-arrival-time copula of memory 1

The IAT 1 is a refined version of SIM: it takes into account additionally the correlation between two successive intervals in the stationary regime by defining a copula. A semi-parametric approach is carried out, wherein the marginal survival functions is assessed in the same way as the SIM method using a Kaplan-Meier estimator, and the copula function is estimated from the Frank copula family, characterized by a single parameter θ asshowninEq.(3). Both negative dependence (θ<0), positive dependence (θ>0) and independence (θ = 0) can be modeled within the Frank family.

C θ (u, v)=- 1 θ log 1+ (e -θu -1)(e -θv -1) e -θ -1 . ( 3
)
The choice of the Frank family is explained as follows. The configurations of the SRP defined in Section 2.3 have been studied. After generating an SRP sequence of length N = 500, 000, several copula families (Gaussian, t, Frank, Gumbel, Clayton and independent copula) are then fitted to the data. It has been found that there exists a certain convergence between the families:

• The degree of freedom ν in the t-copula is systematically estimated to be extremely large, corresponding to the Gaussian copula. Therefore, only the Gaussian copula is maintained in the comparative study

• The parameter of the Gumbel family, θ Gumbel is evaluated to be 1 and that of the Clayton family, θ Clayton is estimated to be 0 + , which corresponds in both cases to the independent copula. As the Gumbel and Clayton copulas do not take into account the negative dependency, both models are not kept for the comparative study.

Since the correlation between the intervals is of interest to the current paper, the independent copula, as well as the Gumbel and Calyton, are excluded from consideration. Gaussian and Frank copula are now the two candidates. In the following, we investigate whether they differ in estimating the tail dependence.

The concept of tail dependence describes the amount of dependence in the lower-leftquadrant tail or the upper-right-quadrant tail of a bivariate distribution. Let F 1 and F 2 be the marginal distributions of variables X 1 and X 2 , respectively. A common measure of tail dependence is given by the so-called upper/lower tail dependence function [START_REF] Frahm | Estimating the tail-dependence coefficient: properties and pitfalls[END_REF]:

λ L (v)=P (F 1 (X 1 ) ≤ v|F 2 (X 2 ) ≤ v)= C(v, v) v , (4)
λ U (v)=P (F 1 (X 1 ) ≥ v|F 2 (X 2 ) ≥ v)= 1 -2v + C(v, v) 1 -v . ( 5
)
An example of the tail dependence functions λ L (v),v ∈ [0, 0.5] and λ U (v),v ∈ [0.5, 1] are shown, respectively, in Figures 4 and5 for the SRP configuration II. The tail dependence of the intervals in SRP is estimated empirically [START_REF] Caillault | Empirical estimation of tail dependence using copulas. Application to Asian markets[END_REF]. It is observed that the Frank (triangle marker) and the Gaussian copulas (circle marker) are both close to the empirical SRP (cross marker). The independent copula is clearly inadequate.

The GOF tests [START_REF] Genest | Goodness-of-fit tests for copulas: a review and a power study[END_REF] are then used to assess to what extent the copula in an SRP resembles a Frank copula or a Gaussian copula. The copula under null hypothesis is Frank or Gaussian. Let α be the level of significance. With infinite data, an SRP copula can never be approximated by a Frank/Gaussian copula, meaning that the null hypothesis will be 100 percent rejected (p-value less or equal to α). When the amount of observations are limited, however, the SRP copula can be estimated very satisfyingly by a Frank/Gaussian copula.

Let α =0.05. An SRP sequence of length N = 20 is generated 1, 000 times, and the corresponding copula (which describes the dependence structure between two successive intervals in an SRP) is tested against Frank/Gaussian, using the Cramer-von Mises statistics [START_REF] Genest | Goodness-of-fit tests for copulas: a review and a power study[END_REF]. The total proportion of an SRP copula being classified as Frank/Gaussian, is recorded in Table 2, column "GOF Frank" and "GOF Gaussian." It is found that the proportion of SRP copula, which "look like" a Gaussian copula, is slightly higher than that for a Frank copula. To sum up, both Frank and Gaussian copulas have similar performance in estimating the tail dependence. In terms of the GOF test, Gaussian is slightly better than Frank. Nevertheless, Frank copula admits an explicit formula which is not possible for Gaussian. This facilitates greatly the inference procedure when missing/censored data is involved, as shown later in the case studies.

Once the marginal functions have been characterized, θ can be estimated by maximum likelihood estimation, and the two predictors can be assessed. In the presence of both CM and PM, the estimation procedure is more complex as it potentially involves doubly censored bivariate data. This issue has been specifically addressed in [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF].

Arithmetic reduction of age with infinite memory

Using imperfect maintenance models such as the ARA ∞ virtual age assumption is an alternate possibility to approximate an SRP. Inference procedures have already been implemented for virtual age models when the system, under CM only, is assumed to be in its stationary regime at the beginning of the observations [START_REF] Dijoux | Statistical inference for imperfect maintenance models with missing data[END_REF] for any given initial intensity. In the following, it is assumed that the initial intensity follows a Weibull distribution, which is a common assumption for aging systems in reliability engineering. For the configurations A and B, three parameters need to be estimated: η and β, respectively, the scale and shape parameter of the initial Weibull intensity and ρ 0 ∈]0, 1], which characterizes the maintenance efficiency. In the presence of CM and PM such as in configuration C, it is also assumed that the PM efficiency follows an ARA ∞ assumption with a parameter ρ 1 , with traditionally ρ 1 >ρ 0 . The parameters of the models can be estimated by maximum likelihood estimation in a similar way as in [START_REF] Dijoux | Statistical inference for imperfect maintenance models with missing data[END_REF]. Let us denote λ(t)=(β/η)(t/η) β-1 the initial Weibull intensity and Λ(t)=(t/η) β the cumulative intensity. First, the likelihood function L a0 associated with the observation of the n first maintenance times and types (X, δ) 1...n and an initial age a 0 can be expressed as in Eq. ( 6). The effective ages in Eq. ( 6) can be obtained by induction given a 0 and using

a i =(1-ρ δi )(a i-1 + x i ). L a0 ((X, δ) 1...n )= n i=1 (λ (a i-1 + x i )) 1-δi × exp n i=1 -Λ(a i-1 + x i )+Λ(a i-1 ) . (6)
Second, as the SRP is assumed to be in its stationary regime at the beginning of the observations, the same assumption is proposed for its approximating model. It implies that a 0 is the realization of the limiting age distribution A ARA ∞ with pdf f ARA A∞ . The resulting likelihood function L associated with the observation of the n first maintenance times and types (X, δ) 1...n and an initial mixing distribution A ARA ∞ is presented in Eq. [START_REF] Cox | On the superposition of renewal processes[END_REF]. The mixing distribution A ARA ∞ has been characterized theoretically in Nguyen et al. [START_REF] Nguyen | Analytical properties of an imperfect repair model and application in preventive maintenance scheduling[END_REF] for CM only and can be derived empirically from intensive simulations considering CM and PM.

L((X, δ) 1...n )= ∞ 0 L a ((X, δ) 1...n )f ARA A∞ (a)da. (7)
Finally, once the parameters of the model estimated, E[X ∞ ] and Corr can be computed numerically, as shown in Nguyen et al. [START_REF] Nguyen | Analytical properties of an imperfect repair model and application in preventive maintenance scheduling[END_REF] and Liu et al. [START_REF] Liu | Approximation of superimposed renewal processes by virtual age models and copulas[END_REF], respectively.

Brown-Proschan model

The BP virtual age assumption is another candidate to approximate an SRP as both converge to a stationary regime. Consider an initial Weibull intensity. Estimation methods have been presented in Lim [START_REF] Lim | Estimating system reliability with fully masked data under Brown-Proschan imperfect repair model[END_REF], Doyen [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF] and specifically in the case of initial steady-state in Dijoux et al. [START_REF] Dijoux | Statistical inference for imperfect maintenance models with missing data[END_REF] for CM only. With only CM, the parameters to estimate are the scale and shape parameters of the Weibull distribution and the CM efficiency p 0 ∈]0, 1]. With both CM and PM, the PM efficiency is assumed to follow a BP assumption with a parameter p 1 . The likelihood function associated with the observation of n maintenance times and types can be derived in three steps. First, the likelihood function L NHPP a0 associated with the observation of a system with initial age a 0 under minimal repair are reminded in Eq. ( 8) using the classic results of non-homogeneous Poisson processes (NHPP) [START_REF] Rigdon | The power law process: a model for the reliability of repairable systems[END_REF].

L NHPP a0 ((X, δ) 1...n )= n i=1 (λ (a 0 + t i )) 1-δi × exp (-Λ(a 0 + t n )+Λ(a 0 )) . (8)
Second, the likelihood function L BP a0 associated with the observation of a system with initial age a 0 and BP efficiencies are derived recursively in Eq. ( 9) using a similar approach as in [START_REF] Dijoux | Statistical inference for imperfect maintenance models with missing data[END_REF].

L BP a0 ((X, δ) 1...n )= n-1 k=1 (1 -p δ k ) ×L NHPP a0 ((X, δ) 1...n + n-1 i=1 L BP a0 ((X, δ) 1...i )p δi n-1 k=i+1 (1 -p δ k ) ×L NHPP 0 ((X, δ) i+1...n ). (9)
Finally, the initial age a 0 is traditionally not known and can be assumed to be the realization of a random variable A BP ∞ , limiting distribution of the effective age in the stationary regime. The resulting likelihood function L BP associated with the observation of the n first maintenance times and types (X, δ) 1...n and an initial mixing distribution A BP ∞ is presented in Eq. [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF].

L BP ((X, δ) 1...n )= ∞ 0 L BP a ((X, δ) 1...n )f A BP ∞ (a)da. ( 10
)
The distribution of A BP ∞ has been expressed theoretically in [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF] with CM only. With CM and PM, the pdf of A BP ∞ can be obtained numerically. Once the parameters of the models are estimated, the predictors can be obtained from a plugged-in version of their theoretical expressions or numerically.

Results

The MSE of the two predictors E[X ∞ ] and Corr are, respectively, plotted in Figures 6 and7 for configuration A, Figures 8 and9 for configuration B and Figures 10 and11 for configuration C.

When no PM is involved (configurations A and B), HPP has the smallest error in the estimation of E[X ∞ ]( s e eFigures 6 and 8). Virtual age models (ARA ∞ and BP) perform better than IAT 1 and SIM. With more than 50 data, the differences between the MSEs given by the tested models are less important.

As for the error of the predictor Corr (see Figures 7 and9), HPP and SIM have the largest error since their intervals are independent. BP is an auto-correlated process with a weak dependence between intervals, but its estimation of correlation is far from accurate. IAT 1 has the smallest MSE when the SRP are homogeneous, which may support the practice of using a Frank copula to approximate the dependence structure in such an SRP. If the SRP differ strongly from each other as in configuration B, ARA ∞ outperforms other models in estimating the correlation between adjacent intervals.

When periodic PM are implemented (Figures 10 and11), the SIM and IAT 1 have a considerable error in E[X ∞ ]. HPP performs the best when estimating E[X ∞ ], but with more data, the advantage of HPP over BP/ARA ∞ is less significant. As for the correlation estimator, the performance of the tested models depend on the data length: with less than 50 data, HPP outperforms the others; otherwise, BP and ARA ∞ are the best. It should be emphasized that the implementation of periodic PM results in a non-null correlation in HPP/SIM sequences.

A BANE NOR CASE STUDY

General Overview and Data Presentation

Signaling is essentially a sophisticated traffic light system for the railway. The complexities of moving trains around such a large network, keeping them safely apart, and allowing for their long stopping distances, means that the signaling system is very complicated and comprises a great many parts. The signals themselves are the line-side pieces of equipment that tell train drivers when it is safe to proceed and what route their train will take. A light signal comprises (1) signal head with light sources, background screen and shadow screens; (2) mast with platform/ladder when needed and (3) devices for controlling the signal with interface to interlocking equipment.

In this section, we focus on the sub-system of the light sources in a signal head, formed by two to five lights of different colors (red, green, yellow and white). One lamp can be either LED or incandescent. A priori, there is no manifest heterogeneity in the lights of the same type, but the maintenance strategy depends on the importance of the lamp: since the failure of a red light which gives the signal "stop" is often more severe than a failure of a green one, preventive maintenances have been performed only on the signals that consists of red light bulbs. Both corrective and preventive maintenances consist in replacing the light bulb, burned or still working, by a new one. Some preventive maintenances, like periodic inspection or cleaning, are planned for other parts of the signal (cables, covering glass, etc.), and are not considered here.

However, the information regarding the maintenance actions is incomplete. When a burned light bulb is replaced, its position in the system is usually not recorded. Thus, the failure history of an individual light bulb is not available and the failures of the system must be treated as a whole. An SRP can, therefore, be used to describe the successive failures of the light sources in a signal.

In the following, three case studies are presented where the proposed five approximation models are fitted to the failures times of dwarf signals. The first study is on a single asset, showing how to use these models to compute the expectation of RUL and evaluating the virtual age; the second one focuses on signals having numerous failure records (more than 7), with the assumption that despite the geographical heterogeneity, these signals could be described by models of the same parameters; the third one studies all the dwarf signals composed by four incandescent light without redundancy (single filament lamp).

The maintenance records provided by Bane NOR contains not only failure histories (CM) but also some invalid PM plans. The earliest date that appeared in a PM record is considered as the beginning of the individual observation window. For the first case studies, the observations start from the first recorded CM; in the second and third case study, the beginning of observation is the earliest date of the planned PM, determined individually for each asset; and the end of observation is the last date of CM, which is March 28, 2019. Thus, for each asset, there is at least one right-censored data (Table 3).

Case Study 1: Asset 012110

The investigated signal is located in Dovre Line, a main national connection between Eastern Norway and Trøndelag and further north for passenger and freight traffic. The consecutive lifetimes are shown in Table 4.

Parameters of fitted models and some reliability indicators are gathered in Table 5.The second column A 0 represents the expected virtual age at the beginning of the observation and needs to be calculated only for imperfect maintenance models, namely ARA ∞ and BP. E[X] is the unconditional expected lifetime in each model. For the models which are not renewal processes, the expected value of the 9th lifetime is smaller than E[X] since the last observed lifetime is relatively large.

The calculation of conditional E[X] under ARA ∞ and BP is recalled here. In ARA ∞ , A 0 is derived along with the parameters, then the consecutive effective virtual ages are 1) and plotted in Figure 12. Denote by A end the virtual age after last observation: the conditional reliability of the next interval is R ARA∞ (t|A end )= exp(-α(t + A end ) β + αA β end ). As for BP, one should first retrieve the maintenance effects, which is described in [START_REF] Lim | Estimating system reliability with fully masked data under Brown-Proschan imperfect repair model[END_REF]. The repair effect of the 1st-7th maintenance actions can be represented by a binary vector B =[0, 0, 1, 0, 1, 0, 1], where 0 signifies a minimal repair and 1 a perfect one. Remark that the maintenance effect after the last CM is non-identifiable in the absence of further data. A 0 is calculated [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF] and the consecutive virtual ages are plotted in Figure 13. After the 8th repair, the system is assumed to be as good as new. The repair after the last observed lifetime has the probability p 0 to be perfect and 1 -p 0 to be minimal. Thus, the reliability of the 9th interval can be expressed as R BP (t|B)=p • exp(-αt β )+(1-p) • exp(-α(t + X end ) β + αX β end ), where X end is the last observed lifetime. It can be observed from Figures 12 and13 that the asset seems to be older under ARA ∞ assumption than it is under the BP model. This is because (1) a low maintenance efficiency was estimated for ARA ∞ , which limits the reduction of age at each repair and (2) three maintenances are estimated to be perfect in the BP model, making the asset relatively young. Consequently, the remaining lifetime predicted by ARA ∞ is smaller than that evaluated by BP (Figure 14).

Case Study 2: the Most Frequently Failed Assets

We investigate the inference and the hazard rate in the 16 most frequently failed assets. One may assume that there exist some clusters of assets in which all members are considered similar and share the same parameters of a certain model. The clustering can be based either on physical/geometrical features, that is, assets located at the same station and composed by lamps of the same type, or on lifetimes features, that is, a bunch of signals whose MTBF are significantly different from the others. The inference procedure for clustered assets is to maximize the sum of the log-likelihood function of all group members. The data set is therefore highly censored. Ignoring the heterogeneity and assume that the failure histories share the same model parameters, the inference results are gathered in Table 8.

Since the Kaplan-Meier estimator is not defined beyond the largest observation which is right-censored, the expected lifetime in the SIM model is calculated by changing the largest observation to a "death." The survival functions given by each model are shown in Figure 16. It can be observed that ARA ∞ and BP differ in terms of the tail distribution (see Figure 17), resulting in a significant difference in the mean lifetime.

The fact that the empirical survival curve crosses its exponential fit once only, and from below, suggests strongly the possibility that the failure distribution has a decreasing failure rate (DFR) [START_REF] Proschan | Theoretical explanation of observed decreasing failure rate[END_REF]. Both ARA ∞ and BP exhibits the DFR via a Weibull shape parameter smaller than 1, and their unconditional survival function is much closer to the Kaplan-Meier estimate than the exponential fit. It has been pointed out in [START_REF] Cha | Some notes on unobserved parameters (frailties) in reliability modeling[END_REF], however, that the observed DFR is possibly related to a heterogeneous population. In our case, the "observed" heterogeneity results directly from some measurable covariates, that is, working environment (average temperature, humidity and precipitation) and usage (total lighted time of the lamps), whereas the "unobserved" heterogeneity includes the variation in the quality or robustness of the lamps. For instance, some of the light bulbs are so robust that they can survive the fluctuation of voltage that could kill the others. Without external damage, their lifetime can be up to several decades.

CONCLUSION

This paper introduces five approximating models of an SRP: two imperfect repair models (ARA ∞ and BP), two renewal models (SIM and HPP) and IAT 1 , constructed by the Kaplan-Meier estimated marginal distribution and a frank copula which captures the dependence structure between successive intervals.

The performances of these models are evaluated by investigating the amplitude of errors of mean interval length and correlations when the above-mentioned approximations are used. Further, their capabilities in the prognosis of RUL is examined using the level-set approach and the Gini index. It is hard to say which model is overall the best: the performances depend on the aging rate as well as the available data amount.

Then, three case studies on signal failure data provided by Bane NOR are conducted to complete the Monte Carlo simulations. It is observed that there may exist significant heterogeneity between the signals: a few assets exhibit the IFR, but when all assets are merged together, the DFR is revealed. The results may justify the current maintenance plan where no PM is implemented. Several limits of the case studies should be emphasized: first, the failure history of any individual signal is short; second, the "burn-in" phase of the lamps has been observed, but this can hardly be modeled by the proposed approximations; third, the DFR nature is related to the heterogeneity in the population, caused by different working conditions or distinctive frequency of usage.

Future research involves a thorough study of the heterogeneity, which plays an important role in survival analysis, especially the failure rate estimation. On the one hand, our Monte Carlo simulations suggest that the when heterogeneous virtual age processes (i.e., ARA ∞ processes that share same aging parameter β and repair effectiveness ρ but differ in the scale parameter α) are merged, traditional ML estimation will under-estimate the aging parameter and over-estimate the repair efficiency; on the other hand, the effect of numerous covariate (environmental conditions, frequency of usage and manufacture) could be analyzed through, for instance, the Cox proportional hazards model or random effect survival models. These studies shall reinforce our understanding about the aging mechanism of signals and are helpful for the optimization of maintenance planning in signaling system. Using the mean value theorem, the right side is reformulated as

λ b s (y) ∞ x R b u (y) - d du Ra(u) du, s ∈ [x, ∞].
When the system is IFR, λ b x (y) ≤ λ b s (y), leading to A.2.

A.3 SRP

We introduce the following corollary when the SRP are identical: (A.

3)

The proof is straightforward. Inequality (A.3) is easily satisfied in the case IFR since λ(t)

∞ t e -Λ(u) du ≤ ∞ t λ(u)e -Λ(u) du = R(t).

Figure 1 .

 1 Figure 1. Illustrative trajectory of a CM process for a dwarf signal.

  Procedure to Measure the Performance of the Models L e tu sc o n s i d e rar i s kα and a trajectory of SRP observations where the parameters of the model and the age of each component are available. At time T i , it is straightforward to compute the confidence interval of smallest amplitude I =[T i + a, T i + b] with confidence level 100(1 -α)% from the actual model. Next, an interval J

Figure 2 .

 2 Figure 2. Level sets in the fast aging rate cases: homogeneous components (configuration II).

Figure 3 .

 3 Figure 3. Level sets in the fast aging rate cases: heterogeneous components (configuration IV).

Figure 4 .

 4 Figure 4. Lower tail dependence function.

Figure 5 .

 5 Figure 5. Upper tail dependence function.

Figure 6 .

 6 Figure 6. Configuration A: MSE of E[X ∞ ].

Figure 7 .

 7 Figure 7. Configuration A: MSE of Corr.

Figure 8 .

 8 Figure 8. Configuration B: MSE of E[X ∞ ].

Figure 9 .

 9 Figure 9. Configuration B: MSE of Corr.

Figure 10 .

 10 Figure 10. Configuration C: MSE of E[X ∞ ].

Figure 11 .

 11 Figure 11. Configuration C: MSE of Corr.

Figure 12 .

 12 Figure 12. Virtual ages in the approximating ARA ∞ model.

Figure 13 .

 13 Figure 13. Virtual ages in the approximating BP model.

Figure 14 .

 14 Figure 14. Survival functions of the inter-failure times of asset 012110.

Figure 15 .

 15 Figure 15. Survival function of the inter-failure times derived from the 16 most frequently failed assets.

Figure 16 .

 16 Figure 16. Survival functions derived from the failure times of all dwarf signals.

Figure 17 .

 17 Figure 17. Tails of the survival function under ARA ∞ and BP models.

X

  . Liu et al. Note Ra(t)=P (Xn ≥ t|A n-1 = a)a n dR b u (t)=P (X n+1 ≥ t|An =(1-ρ)(a + u)). The derivative of R(x, y) with respect to x, y and mixed derivative are, respectively: The derivative of survival function R b x (t)is the product of failure rate and survival: -(d/dy)R b x (y)=λ b x (y)R b x (y). Therefore, R is RR 2 if and only if

Corollary A. 1 :

 1 Note R, f the survival function and density function of an interval in the superposed RP. Define φ(t)= ∞ t R(u)du. Then, two successive intervals in the SRP in stationary regime is RR 2 -dependent if and only ifφ(t)f (t) ≤ R(t) 2 .

Table 1 .

 1 Gini coefficients(10 -3 )

		HPP	SIM	IAT 1	ARA∞	BP
	I	9.1	8.9	8.5	8.2	8.6
	II	119.2	119.5	137.0	89.4	144.0
	III	6.5	7.2	6.2	3.0	5.9
	IV	222.7	230.1	116.1	133.8	204.5
	V	2.8	4.3	4.3	3.0	5.6

Algorithm 1

 1 Empirical MSE computation 1: Select a SRP configuration 2: Compute the predictors e=E[X ∞ ] and r=Corr 3: for N ∈{10, 20, 50, 100, 200} //N=length of a trajectory

	4:	for k=1 : 5000 //kth sample
	5:	Simulate a SRP trajectory Hist of length N
	6:

Table 2 .

 2 Estimated copula parameters and results of test of GOF

	SRP	θ in Frank copula	ρ in Gaussian copula	GOF Frank	GOF Gaussian
	I	-0.7372	-0.1234	0.949	0.966
	II	-2.3652	-0.3634	0.926	0.952
	III	-0.6137	-0.1046	0.955	0.963
	IV	-2.0441	-0.3151	0.926	0.946

Table 3 .

 3 A summary of the case studies

		Object	Inference	Censoring	No. asset
	Case study 1	Asset 012110	Individual	No	1
	Case study 2	Signals having 7 + failures	Grouped	Yes	16
	Case study 3	All dwarf signals	Grouped	Yes	1608

Table 4 .

 4 Inter-failure times of asset 012110

	133	10	283	763	19	378	203	920

Table 5 .

 5 Parameters of the models fitted to asset 012110

		Parameters	A 0	E[X]	Expected value of 9th lifetime
	ARA∞	α =3.34 × 10 -10
		β =2.9737	1717.1	321.38	239.83
		ρ =0.1577
	BP	α =3.27 × 10 -10
		β =3.3571	460.77	321.91	293.36
		p =0.4113
	HPP	339.08	0	339.08	339.08
	IAT 1	θ = -1.32	0	281.52	210.49
	SIM		0	281.52	281.52

Table 6 .

 6 Reliability indicators for the 16 most frequently failed assets

		Parameters	A 0	E[X]
	ARA∞	α =2.16 × 10 -12
		β =3.5363	1706.1	502.97
		ρ =0.2277
	BP	α =2.55 × 10 -7
		β =2.2417	575.07	509.63
		p =0.4698
	HPP	460.27	0	460.27
	IAT 1	θ =0.1365	0	440.12
	SIM		0	440.12

Table 7 .

 7 Distribution of the number of recorded CM

	N o . . f a i l u r e0	1	2	34567	8	9	1 0	1 1
	No. asset	821	419	178	95	51	17	11	8	2	4	0	2
	4.4. Case Study 3: All Dwarf Signals
	1,608 assets are investigated, half of which has no CM record. Table 7 gathers the
	distribution of number of CM of these assets.

Table 8 .

 8 Parameters and reliability indicators for all 1,608 dwarf signals

		Parameters	A 0	E[X]
	ARA∞	α =0.2619
		β =0.1948	13,448	474,963 (1,300 years)
		ρ =0.9723
	BP	α =0.0043
		β =0.5625	0	26195 (72 years)
		p =1.0000
	HPP	6227.1	0	6227.1 (17 years)
	SIM		0	2576.1 (7 years)

Acknowledgments

The authors thank the anonymous reviewers for their valuable comments. The Ph.D. project of Xingheng Liu is funded by the Ministère de L'enseignement Supérieur et de la Recherche (French Minister of Education and Research) and the European Regional Development Fund, through contract with the Universitéd e Technologie de Troyes.

A. PROOF OF PROPOSITION 2.1 A.1 ABOUT RR 2 FUNCTION

A pair of real-valued random variables (X 1 ,X 2 ) and its density function f (•, •) are called reverse rule of order 2 (RR 2 )( [START_REF] Lehmann | Some concepts of dependence[END_REF]) if

whenever x 1 >x 2 and y 1 >y 2 . According to Karlin and Rinott [START_REF] Karlin | Classes of orderings of measures and related correlation inequalities II. Multivariate reverse rule distributions[END_REF], this is equivalent to

The joint survival of two successive intervals in ARA in the stationary regime, conditioned on the previous virtual age satisfies Eq. (A.2). The proof is given below:


Jad Botros
email: jad.botros@utt.fr

Farah Mourad-Chehade
email: farah.chehade@utt.fr

David Laplanche
email: david.laplanche@hcs-sante.fr

Detection of Heart Failure Using a Convolutional Neural Network via ECG Signals

Keywords: Deep Learning, Binary Classification, Convolutional Neural Network, Heart Failure, Electrocardiogram

Heart failure (HF) is a chronic heart condition that increases mortality, morbidity, and healthcare costs. The electrocardiogram (ECG) is a noninvasive and straightforward diagnostic tool that can reveal detectable changes in HF. Because of their small amplitude and duration, these changes can be subtle and potentially misclassified during manual interpretation or when analyzed by clinicians. This paper reports a 7-layer deep convolutional neural network (CNN) model for HF automatic detection. The proposed CNN model requires only minimal preprocessing of ECG signals and does not require any engineered features. The model is trained and tested using an unbalanced and a balanced datasets extracted from the MIT-BIH and the BIDMC databases, achieving an accuracy of 99.73%, a sensitivity of 99.58%, and a specificity of 99.83% when the dataset is unbalanced and an accuracy of 99.26%, a sensitivity of 99.37%, and a specificity of 99.12% when the dataset is balanced.

I. INTRODUCTION

Heart failure (HF) is the inability of the heart to pump enough blood to meet metabolic demands or to supply venous return. It is a clinical syndrome caused by muscular tissue injuries in the heart. Ischemic heart disease, hypertension, and diabetes are the leading causes of these injuries. As the heart fails, patients experience a variety of symptoms such as dyspnea due to pulmonary congestion, edema and ascites due to impaired venous return. Common constitutional symptoms include nausea, loss of appetite, and fatigue. Several compensatory mechanisms occur as the failing heart attempts to maintain adequate function, such as increasing cardiac output, increasing ventricular volume and wall thickness, and maintaining tissue perfusion with augmented mean arterial pressure, among other things. Although these mechanisms may be beneficial in the early stages of HF, they all eventually lead to a vicious cycle of worsening HF [START_REF] Kemp | The pathophysiology of heart failure[END_REF].

According to statistics, HF affects 2 to 3% of the population in developed countries: approximately 15 million Europeans are affected by this syndrome [START_REF] Braunschweig | What are the costs of heart failure?[END_REF]. Globally, more than 26 million people have been diagnosed with HF, creating a major public health issue as well as a significant economic burden [START_REF] Ponikowski | Others Heart failure: preventing disease and death worldwide[END_REF]. The prevalence of this syndrome has been increasing in recent years and is expected to rise further as the population ages and even patients with cardiac conditions, which may trigger the development of HF, live longer. Despite the advances in HF care, the disease still has a high morbidity and mortality rate: nearly half of all patients diagnosed with HF die within five years. HF, on the other hand, is associated with frequent hospitalization: in Europe, HF accounts for 5% of all acute hospital admissions [START_REF] Braunschweig | What are the costs of heart failure?[END_REF]. HF care accounts for 1 to 2% of total health-care costs in developed countries. The cost of hospitalization accounts for the majority of total spending [START_REF] Berry | Economics of chronic heart failure[END_REF].

A combination of signs and symptoms, corroborated by tests, is required for the diagnosis of HF. The electrocardiogram (ECG) is a non-invasive measurement used to record the electrical activity of the heart. HF has an effect on such a signal. However, visual evaluation of the ECG signal obtained from the patient takes time, and manual interpretation is subject to inter-observer variability. The changes that hit the ECG signal are of a differential nature, making them difficult to be detected with the naked eye. This highlights the importance of using Computer Aided Diagnosis (CAD) systems to detect such abnormalities, because of their improved accuracy and reliability. A typical CAD system consists of four key procedures: signal pre-processing, extraction of unique features, selection of significant features, and classification. This paper proposes a novel model based on deep learning to diagnose HF with a short-time record of ECG signals.

The proposed model is a deep Convolutional Neural Network (CNN) having seven layers, with minimal pre-processing requirements. The main advantage of using a CNN over traditional machine learning methods is that it automatically detects important features without any human intervention. With its seven layers, the proposed architecture is smaller than currently available architectures in the literature, and thus it is less consuming in memory and computation time. This approach is illustrated using real ECG signals acquired from ambulatory holter and ECG recorders databases. Evaluation metrics demonstrate better performances, even with 2s ECG segments as input.

The rest of the paper is organized as follows. Section II presents related work. Section III describes the database and the proposed method. Section IV presents and discusses the obtained results. Finally, Section V provides the conclusion.

II. RELATED WORK

Many machine learning and deep learning techniques have been proposed in the literature to address the HF diagnosis issue. Asyali investigated the discrimination power of nine commonly used long term Heart Rate Variability (HRV) measures in his study [START_REF] Asyali | Discrimination power of long-term heart rate variability measures[END_REF]. He then used a Bayesian classifier to classify ECG signals as Normal or HF based on the standard deviation of NN intervals or SDNN. The obtained accuracy equals 93.24%, with sensitivity and specificity of 81.82% and 98.1% respectively. Melillo et al. distinguished between Normal and HF patients using a combination of HRV measures fed into a Classification and Regression Trees model (CART) [START_REF] Melillo | Discrimination power of long-term heart rate variability measures for chronic heart failure detection[END_REF]. The classification accuracy, sensitivity and specificity were 96.36%, 89.74% and 100% respectively. Liu et al. described a method for detecting HF using three other combinations of HRV measures and a Support Vector Machine (SVM) classifier [START_REF] Liu | A new approach to detect congestive heart failure using short-term heart rate variability measures[END_REF]. The model was perfect in terms of accuracy, sensitivity and specificity. Masetic et al. used feature extraction and classification to achieve the same goal [START_REF] Masetic | Congestive heart failure detection using random forest classifier[END_REF]. The autoregressive burg method was used to extract the features, which were then fed into a Random Forest classifier, which achieved 100% accuracy, sensitivity, and specificity. Data pre-processing is an essential step in machine learning techniques because the quality of the data, and thus the information that can be extracted from it, directly affects the network's learning ability. Another important step in machine learning is feature selection. The selection process can be time-consuming and labor-intensive. Deep learning models have been used to optimize the performance of a CAD system in order to avoid the pitfalls of traditional machine learning. Chen et al. used RR interval segments to detect HF using a deep learning model based on sparse auto-encoders (SAE) [START_REF] Chen | A CHF detection method based on deep learning with RR intervals[END_REF]. The SAE were used to extract some features, which were then fed to a Deep Neural Network (DNN) for classification. The results showed a 72.44% accuracy, a 50.93% sensitivity, and a 80.93% specificity. Wang et al. presented a deep learning method based on Long Short Term Memory (LSTM) to classify ECG signals as Normal or HF [START_REF] Wang | CHF Detection with LSTM Neural Network[END_REF]. A 500-point RR intervals segment was fed into the network. The model's accuracy was 84.91%, its sensitivity was 75.49%, and its specificity was 90.06%. Acharya et al. proposed an 11-layer Convolutional Neural Network (CNN) for categorizing ECG signals as Normal or HF [START_REF] Acharya | Deep convolutional neural network for the automated diagnosis of congestive heart failure using ECG signals[END_REF]. The accuracy was 98.97%, the sensitivity was 98.87%, and the specificity was 99.01%. Wang et al. used LSTM and an inception module to detect HF. For an input of 500 points RR interval segment, the network achieved 86.42% accuracy, 74.91% sensitivity, and 91.21% specificity [START_REF] Wang | Detection of congestive heart failure based on LSTM-based deep network via short-term RR intervals[END_REF]. Padmavathi et al. discussed an 11-layer CNN approach for detecting

III. MATERIALS AND METHOD

A. Database Description

The ECG signals used in this retrospective study are obtained from two publicly accessible databases: the MIT-BIH normal sinus rhythm (MIT-BIH NSR) database [START_REF] Goldberger | Phys-ioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals[END_REF] and the BIDMC congestive heart failure (BIDMC CHF) database [START_REF] Baim | Survival of patients with severe congestive heart failure treated with oral milrinone[END_REF]. The former database contains eighteen long-term ECG signals from five healthy men aged 26 to 45 and thirteen healthy women aged 20 to 50. The ECG signals are sampled at 128 Hz and cover 20h of recordings. The latter database contains fifteen long-term ECG signals from eleven men aged 22 to 71 and four women aged 54 to 63 who are diagnosed with heart failure. The HF ECG signals are sampled at a frequency of 250 Hz, also with 20h of recordings. To ensure that all ECG signals are sampled at the same frequency, the MIT-BIH NSR recordings are upsampled to 250 Hz. The signals are then segmented into 2s fragments. This leads to 36 000 2-s ECG segments per signal. Following that, the segments are regularized using Z-score normalization to ensure that they all contribute equally to the analysis. Table II The generated dataset lacks balance, with unequal representation of the two classes (60 % normal and 40 % HF). A balanced dataset is created to avoid the problem of data imbalance by taking equal numbers of ECG segments from the MIT-BIH NSR and the BIDMC CHF databases. The resulting dataset is composed of 240 000 2s-ECG fragments, with 120 000 for each class, normal or HF. More fragments could be taken, at the cost of higher computation times and memory consumption. In order to involve equally all subjects for each class, with 18 normal subjects and 15 HF patients: out of the 18 normal subjects, 15 are randomly chosen, and 8000 2-s ECG samples are randomly selected per patient. The obtained dataset will be used in the following for training and testing the model.

B. Method

This paper proposes a 7-layer CNN. The main advantage of using a CNN over its predecessors is that it can detect important features automatically without human intervention. Its low reliance on pre-processing reduces the need for human effort while expanding its capabilities. The model consists of two convolutional layers, two max-pooling layers, and three fully-connected layers. The convolutional layer is made up of a fixed number of filters and is used to extract the maps of the features. Convolutional layers perform the convolution of different filters with the input signal according to equation 1 in which Y is the result of the convolution of the signal x by the kernel f , n and k are positional indexes, and N is the size of the signal. The max-pooling layer, which takes the maximum value in a specific filter region, is used to reduce the network's dimensionality. The fully connected layer, which is usually placed before the output layer, is used to aggregate data from the final feature map and generate the final classification. Fig. 1 shows the general structure of the proposed network.

Y n = N -1 k=0 x k f n-k , (1)
The first layer takes the 2-s ECG segment of 500 points as input. It is a convolutional layer with five 1 × 13 filters applied with a stride of 1 i.e. five feature maps are generated by convolving the different filters with the input 500 points ECG signal. The second layer is a max-pooling layer with a pool size of two and a stride of 4. This layer reduces the dimensions of the features maps by convolving a 1 × 2 filter with each of the previously generated feature maps. As a result, the number of parameters to learn is reduced, as is the amount of processing in the network. Consequently, the model is more robust to changes in feature position in the input. Following that, another convolutional layer of ten filters of size 1 × 9 each is applied with a stride of 1. This set of filters is applied to the dimensionally reduced feature maps extracting a higher-level features. The fourth layer is The model weights are initialized using Glorot uniform initialization, and they are updated using backpropagation with a batch size of 10. The model is built up over a period of sixty epochs. Let y i be the label of segment i of the training databaset, that is y i = 1 if the segment i is taken from a HF signal and y i = 0 otherwise. Consider pi the predicted probability that the segment i is HF, obtained at the output of the network. For the binary classification problem, the binary cross entropy function is used to calculate the model's loss, as shown in the following equation:

L(q) = -1 M Σ M i=1 y i . log(p i ) + (1 -y i ). log(1 -pi ),
where M is the total number of segments taken in the training phase and q is the index of the epoch. The cross entropy computes a score that represents the mean difference between the actual and the predicted values. The score is to be minimized, where the 0 value is a perfect cross entropy.

IV. RESULTS AND DISCUSSION

The model is validated using stratified 10-fold cross validation over the dataset of 208 000 segments. The dataset is divided into ten equal parts with an equal percentage of each class in each part, and the training and the validation are performed ten times, each time using 9 folds for the training and the remaining one for the validation. Following this stratified 10-fold cross-validation procedure, a blindfold validation of normal and HF ECG segments from 4 patients (32 000 segments) with complete holdout is conducted in order to further assess the proposed model. This process is repeated 10 times and Tables III and IV display and summarize the results from this fully held-out test set considering the unbalanced and balanced datasets respectively.

The test set of fully unseen patients shows that the proposed classifier achieves an overall average classification accuracy of 99.263%, correctly categorizing 99.374% of normal ECG fragments and 99.124% of ECG fragments exhibiting indications of HF. Three evaluation metrics are used in this paper to assess the performance of the proposed model that are accuracy, sensitivity and specificity, as listed below:

Accuracy = T P + T N T P + F N + T N + F P Sensitivity = T P T P + F N Specif icity = T N T N + F P
where T P , T N , F P and F N denote the true positive, the true negative, the false positive and the false negative rates.

Tables III and IV summarize the aforementioned measures for the executed ten folds. As can be seen, the presented model achieves an average accuracy of 99.73%, with an average sensitivity of 99.58% and an average specificity of 99.83% with an unbalanced dataset and an average accuracy of 99.26%, with an average sensitivity of 99.37% and an average specificity of 99.12%, with respectively low standard deviations with a balanced dataset. The number of epochs chosen is sufficient for the loss function to converge and reach a stable minimum. The model's average loss is 0.024 when the dataset is unbalanced and 0.08 when the dataset is balanced. The code was written in Python and implemented in Google Colab.

In order to illustrate the performances of the proposed method, a comparison is established between the achieved results and those of other HF diagnosis methods using the same databases. Three existing methods are thus considered. The authors of [START_REF] Melillo | Discrimination power of long-term heart rate variability measures for chronic heart failure detection[END_REF] explored the discrimination power of various Heart Rate Variability (HRV) measures and chose three combinations of these measures to form the input to a CART classifier. The method in [START_REF] Acharya | Deep convolutional neural network for the automated diagnosis of congestive heart failure using ECG signals[END_REF] used an 11-layer CNN to classify 2s ECG fragments, with minimal pre-processing, into N or HF. The authors of [START_REF] Padmavathi | Heart Disease Recognition from ECG Signal Using Deep Learning[END_REF] upsampled the ECG signals to a common frequency, filtered them using empirical mode decomposition, and performed R-peak detection prior to beats segmentation. For classification, the beats were fed into an 11-layer CNN. Table V shows the test evaluation metrics of methods of the literature compared to the proposed model. The results of the method in [START_REF] Acharya | Deep convolutional neural network for the automated diagnosis of congestive heart failure using ECG signals[END_REF] given in Table V are obtained with an unbalanced dataset in the higher row. When balanced data are considered, the accuracy drops to 94.40%, as do the sensitivity and the specificity, which drop from 98.87% to 94.68% and from 99.01% to 94.12%, respectively, as shown in the lower row. The results show that the proposed model performs better than existing methods. A higher specificity is achieved with method of [START_REF] Melillo | Discrimination power of long-term heart rate variability measures for chronic heart failure detection[END_REF], however it is at the cost of the sensitivity and the accuracy that are lower than the obtained one with the proposed model. The use of deep learning in the presented method, along with an efficient training and a suitable architecture, ensures high performances compared to traditional machine learning methods. Indeed, deep learning capabilities enable more reliable signal abstraction in highdimensional space without the need for human intervention. When compared to the deep models architectures in [START_REF] Acharya | Deep convolutional neural network for the automated diagnosis of congestive heart failure using ECG signals[END_REF] and [START_REF] Padmavathi | Heart Disease Recognition from ECG Signal Using Deep Learning[END_REF], the proposed model is smaller, with 7 layers instead of 11, and faster when training time is taken into account.

V. CONCLUSION

An automated ECG classification system for the detection of heart failure is proposed in this paper. The proposed model consists of a 7-layer CNN and is fully automatic, with no preprocessing or feature selection stages required. Normal ECG recordings are obtained from the MIT BIH NSR database, and heart failure ECG recordings are obtained from the BIDMC CHF database. The two databases are freely accessible on the internet. To avoid data imbalance issues, the dataset is created with 120000 normal ECG segments and 120000 ECG segments from heart failure patients. The presented CNN classifier, which achieved high metrics, plays an important role in the detection and classification of ECG signals, as evidenced by its accuracy of 99.73%, sensitivity of 99.58% reflecting the true positive rate, and specificity of 99.83% which represents the true negative rate, when the dataset is unbalanced, and its accuracy of 99.26%, sensitivity of 99.37%, and specificity of 99.12%, when the dataset is balanced. Future work may include the fusion of various medical data, such as times series, punctual measurements and clinical notes, to ensure an efficient monitoring of the state of a subject over time and to offer a more precise classification with a stratification of the patients into different levels of heart failure. The problem of classification would then be defined as a multi-class hybrid classification, combining data of different natures.

Fig. 1 .

 1 Fig. 1. General structure of the proposed model

  summarizes the characteristics of each database.

TABLE III TEN

 III FOLDS EVALUATION METRICS (UNBALANCED DATASET)

	Fold	Accuracy (%)	Sensitivity (%)	Specificity (%)
	1	99.70	99.62	99.75
	2	99.94	99.96	99.91
	3	99.98	99.97	100
	4	99.95	99.93	99.96
	5	99.83	99.53	99.80
	6	99.91	99.97	99.82
	7	98.90	97.50	99.89
	8	99.40	99.56	99.27
	9	99.99	99.98	99.99
	10	99.95	100	99.90
	Average (%)	99.73	99.58	99.83
		TABLE IV
	TEN FOLDS EVALUATION METRICS (BALANCED DATASET)
	Fold	Accuracy (%)	Sensitivity (%)	Specificity (%)
	1	96.74	99.95	94.25
	2	99.85	99.68	100
	3	99.91	99.80	100
	4	99.98	99.96	100
	5	99.05	99.27	99.12
	6	99.97	99.98	99.96
	7	98.44	96.22	100
	8	99.74	99.13	100
	9	99.69	99.82	99.58
	10	99.15	99.93	98.33
	Average (%)	99.26	99.37	99.12


Edgar Ricardo

Silva Russi
email: edgar.silva_russi@utt.fr

Nacima Labadie
email: nacima.labadie@utt.fr

Caroline Prodhon
email: caroline.prodhon@utt.fr

E R Silva

The Grey Zone Two-Echelon Vehicle Routing Problem with Customerto-Parcel Locations and Low-Pollution Vehicles for Inner-City Logistics

Keywords: Echelon vehicle routing problem, Synchronization, Customer to parcel, City logistics, Last-mile delivery, Sustainable logistics

This study addresses the two-echelon vehicle routing problem with grey zones and Customer-to-Parcel (C2P) stations. This problem arises in the search for new sustainable delivery schemes for e-commerce and last-mile distribution in urban areas and aims at reducing last-mile transportation costs. This study proposes a literature review on the subject, and a mixed integer linear programming (MILP) formulation to model and solve small instances of the described problem.

Introduction

Last-mile logistics has gained greater importance as the demand for Business-2-Consumer products in urban areas has increased with the strong growth of ecommerce driven by a wider access to internet services, both for customers and for companies that buy and sell their products through this channel [START_REF] Ismail | Synthèse du problème de routage de véhicules[END_REF]. This generates a huge potential for consolidation and coordination of distribution flows, which play a key role in multi-echelon distribution systems [START_REF] Lewczuk | Vehicle routing in an urban area: Environmental and technological determinants[END_REF] and can help to reduce traffic volume by improving the use of transportation resources [START_REF] Gonzalez-Feliu | The two-echelon capacitated vehicle routing problem[END_REF]. To address last-mile particular challenges, the transition from the use of Internal Combustion Engine Vehicles (ICEV) to Alternative Fuel Vehicles (AFV), the implementation of access restrictions into urban areas and the integration of Customer-to-Parcel stations (e.g. lockers and pick-up points) have become commonly adopted measures. Therefore, we introduce the Grey Zone Two-Echelon Vehicle Routing Problem with Customerto-Parcel Locations and Low-pollution Vehicles for Inner-city Logistics. Different from the classical Two-Echelon Vehicle Routing Problem (2E-VRP), clients may be served via either C2P stations or home delivery in a specified time window; furthermore we consider satellites as possible C2P stations and grey-zone customers on city borders. The contributions of this paper are threefold, as described in the following sections: First, a review of the existing literature is given in Section "Literature Review", then Section "Problem Description and Mathematical Formulation", defines the problem and presents a mathematical formulation to model it. Finally, the results and conclusions are given in Sections "Results" and "Conclusion and Future Research".

Literature Review

Regarding innovative notions into 2E-VRP models intended for last-mile logistics, [START_REF] Anderluh | Multi-objective optimization of a two-echelon vehicle routing problem with vehicle synchronization and 'grey zone' customers arising in urban logistics[END_REF] introduced the notion of grey zones in their multi-objective problem denoted 2eVRPSyn. Customers in such zones are accessible by fleets of the two levels, giving more flexibility to their model. They proposed a large neighborhood search embedded in a heuristic rectangle/cuboid splitting to efficiently solve the problem. Moreover, [START_REF] Qiu | The two-echelon production routing problem with cross-docking satellites[END_REF] introduced the 2-Echelon Production Routing Problem with Crossdocking Satellites (2E-PRPCS), it features production and inventory decisions and it is solved by a Branch-and-Cut algorithm coupled with a first-solution matheuristic while [START_REF] Grangier | An adaptive large neighborhood search for the two-echelon multiple-trip vehicle routing problem with satellite synchronization[END_REF] introduced a new variant called the Two-Echelon multiple-trip Vehicle Routing Problem with Satellite Synchronization (2E-MTVRP-SS) with time windows and service times.

On the other hand, to address e-grocery last-mile challenges, [START_REF] Leyerer | Decision support for urban E-grocery operations completed research[END_REF] proposed a Decision Support System (DSS) for a model in which the second echelon deliveries are performed through either C2P stations or electric cargo bikes. [START_REF] Liu | Two-echelon vehicle-routing problem: Optimization of autonomous delivery vehicle-assisted E-grocery distribution[END_REF] developed an state of the art Hybrid Immune Algorithm (HIA) to solve a multi-objective Two-Echelon Location-Routing Problem with Mixed Vehicles and Mixed Satellites (2E-LRP-MVMS) in which the second echelon deliveries are performed through either C2P stations or Autonomous Delivery Robots (ADRs). Further, [START_REF] Liu | Two-echelon vehicle-routing problem: Optimization of autonomous delivery vehicle-assisted E-grocery distribution[END_REF] modeled the Two-Echelon Vehicle Routing Problem with Mixed Vehicles (2E-VRP-MV) where Autonomous Delivery Vehicles (ADVs) serve the clients in the second echelon. To solve it, a two-step clustering-based hybrid Genetic Algorithm and Particle Swarm Optimization (C-GA-PSO) algorithm are presented. Moreover, as environmental factors had gained more importance in last-mile delivery schemes, [START_REF] Anıl | Variable neighborhood search for the two-echelon electric vehicle routing problem with time windows[END_REF] introduced the 2-Echelon Electric Vehicle Routing Problem with Time Windows (2e-EVRP-TW). The model takes into account charging stations, charging times and time windows. To solve it, a heuristic based on Clarke and Wright and a Variable Neighborhood Search (VNS) are presented. For supplementary information on 2E-VRP refer to [START_REF] Sluijk | Two-echelon vehicle routing problems: A literature review[END_REF].

Furthermore, it was until 2018 when [START_REF] Zhou | A multi-depot two-echelon vehicle routing problem with delivery options arising in the last mile distribution[END_REF] introduced the Multi-Depot Two-Echelon Vehicle Routing Problem with Delivery Options for Last Mile Distribution (MD-TEVRP-DO) that the concept of delivery options through C2P stations was incorpo-ratedintoa2E-VRP . Moreover, [START_REF] Redi | A simulated annealing algorithm for solving two-echelon vehicle routing problem with locker facilities[END_REF] proposed a Simulated Annealing Algorithm to minimize the transportation cost in their conception of 2E-VRP with locker facilities; while [START_REF] Enthoven | The two-echelon vehicle routing problem with covering options: City logistics with cargo bikes and parcel lockers[END_REF] introduced the Two-Echelon Vehicle Routing Problem with Coverage Options (2E-VRP-CO) where last mile deliveries are performed via C2P stations and cargo bikes which is then solved thanks to an Adaptive Large Neighborhood Search (ALNS). Similarly, [START_REF] Anderluh | Synchronizing vans and cargo bikes in a city distribution network[END_REF] considered the synchronization of cargo bikes and vans in their two-echelon model solved by a Greedy Randomized Adaptive Search heuristic with path re-linking for real-world data of the city of Vienna. Later, [START_REF] Yu | Design of a two-echelon freight distribution system in last-mile logistics considering covering locations and occasional drivers[END_REF]i n t r oduced the concept of occasional drivers in a Two-Echelon Vehicle Routing Problem with Time Windows, Coverage Options, and Occasional Drivers (a crowd-shipping concept first introduced by [START_REF] Macrina | The vehicle routing problem with occasional drivers and time windows[END_REF] that can provide potential advantages by its implementation). For further information on the new developments in last-mile logistics, see [START_REF] Boysen | Last-mile delivery concepts: A survey from an operational research perspective[END_REF].

Problem Description and Mathematical Formulation

The problem that arises involves two different fleets of vehicles making deliveries in two different echelons. In our particular case, we consider that the vehicles start and end their routes at their respective depot. ICEVs deliver the goods from the depot of the first echelon, where the entire stock of goods is located, to the first level customers within their time windows and/or to the exchange points, also called satellites. These last are also used as Customer-to-Parcel (C2P) stations where customers can pick up their products directly. AFVs start from the second echelon depot without any goods in cargo so they must immediately meet ICEVs at the satellites where they will pick up the goods to start their deliveries. Then the AFVs serve both direct delivery customers within their respective time windows and second-level C2P stations where customers can retrieve their goods. Due to the costs associated with the use of vehicles and satellites, the waiting time of vehicles at these locations is limited to a specific value and is minimized by the inclusion of the waiting costs in the objective function (54.1) of the MILP model. For this reason, the arrival of the ICEVs and AFVs at the satellites must be approximately at the same time so that the transfer of the goods is as fast as possible and both the ICEV and the AFV can continue their deliveries. Once an AFV has finished its deliveries it can rejoin an ICEV at a satellite to retrieve new goods and continue a new route as long as the maximum route time represented by the delivery operator's working day is not exceeded. Likewise, the same satellite can be used by different AFVs as long as the ICEV's maximum waiting time at the satellite is not exceeded. In the traditional formulation of the 2E-VRP the customers of each echelon are pre-allocated to the first or second echelon, however as [START_REF] Anderluh | Multi-objective optimization of a two-echelon vehicle routing problem with vehicle synchronization and 'grey zone' customers arising in urban logistics[END_REF]shows, this pre-allocation can lead to poor quality results because customers close to the satellites could be served by either an ICEV or an AFV. Thus, [START_REF] Anderluh | Multi-objective optimization of a two-echelon vehicle routing problem with vehicle synchronization and 'grey zone' customers arising in urban logistics[END_REF] introduced the grey zone, an area in which customers are not predefined and can therefore be served by any vehicle. In this model, grey zones will be taken into account only for home delivery customers located on the city borders. On the other hand, C2P stations are increasingly used to eliminate the risk of unattended deliveries and to optimize delivery costs thanks to the possibility of accumulating multiple demands at the same station. Similarly, these stations benefit customers by allowing them to pick up their deliveries at times that are convenient for them at stations close to their homes. Each of these stations has a specific capacity and coverage radius. The mathematical model is inspired on [START_REF] Anderluh | Multi-objective optimization of a two-echelon vehicle routing problem with vehicle synchronization and 'grey zone' customers arising in urban logistics[END_REF] but it considers time windows and waiting time on client nodes, C2P stations for both levels, and location decisions for satellites and stations (Table 54.1).

The objective function (54.1) takes into account the fixed cost of vehicles, the cost per distance, the vehicle operating cost per hour for each level, the cost of opening satellites and C2P stations, and the waiting costs at satellites and client nodes. The classic restrictions of the two-echelon models were formulated (i.e. restrictions on flow, exit and return to the depots, use of vehicles, and service to home-delivery customers for each level).

min i∈v all j∈v all k∈F [(T (i, j ) + ST (i )) * CT (k) + (Dis(i, j) * CD(k))] * X (i, j, k)+ i∈V all k∈F W (i, k) * CT (k) + e∈E j∈v all k∈F X (e, j, k) * FC(k) + i∈c2 p Y (i ) * CU (i ) (54.1)
Also, temporal constraints have been considered to guarantee the service to the clients within their time window. Similarly, for node routing (54.2), the time between nodes, the service time, and the waiting time at satellites and at client's nodes (in case of arrival before the lower limit of the time window) have been considered. Finally, (54.3) limits the maximum duration of the routes.

T ( j, k) ≥ (T (i, k) + Time(i, j ) + ST (i) + W (i, k)) -M * (1 -X (i, j, k)) ∀i ∈ v 1,2 , j ∈ v 1,2 , k ∈ F 1,2 (54.2) T (e f , k) ≤ TM AX ∀e f ∈ E f , k ∈ F (54.3)
Subsequently, constraints (54.4)-(54.6) guarantee that the vehicle load contains the demands of both home-delivery customers and C2P stations. Furthermore, capacity constraints have been formulated for all vehicles and C2P stations.

U ( j, k) + D( j ) ≤ U (i, k) + M * (1 -X (i, j, k))∀i ∈ v all , j ∈ chd, k ∈ F (54.4) U ( j, k) + c∈cloc 2 S(c, j ) * D(c) ≤ U (i, k) + M * (1 - s X (i, j, k)) ∀i ∈ v all , j ∈ loc, k ∈ F 2 (54.5)
Q(k)/QLoc(c2p)
Capacity of vehicle (k) and capacity of the (c2p)station

CT/CD(k)

Cost of using the vehicle (k) per unit of time and distance

CF(k)

Fixed costs of using the vehicle (k)

CU/Radius(c2p)

Opening cost and Covering radius of (c2p)station

ST/D(i)

Service time and demand in node (i)

EA/LA(chd)

Earliest and latest arrival time to home delivery clients (chd) Va r i a b l e s

X(i,j,k)

Binary variable, 1 if the arc (i,j)iscrossedbythe vehicle(k)

T(i,k) Arrival time of vehicle (k) at node (i) Y(c2p) Binary variable, 1 if the (c2p) station is opened, 0 otherwise U(i,k) Load of vehicle (k) at node (i) W(i,k) Waiting time of vehicle (k) at node (i) S(cloc,c2p) Binary variable, 1 if the (c2p) station serves the client (cloc) U ( j, l) + U ( j, k) + c∈cloc 1 S(c, j ) * D(c) ≤ U (i, k) + M * (1 -X (i, j, k)) ∀i ∈ v all , j ∈ sat , , k ∈ F 1 , l ∈ F 2 (54.6)
Moreover, (54.7) ensure that if a satellite is used by an ICEV, it must also be used by an AFV. Likewise, the waiting times are calculated for home delivery customers and satellites, whereas (54.8) determine the maximum waiting time allowed on satellites. Finally, (54.9) ensure that vehicles do not contain cargo when returning to their respective depot at the end of their routes, and (54.10) guarantee that AFVs do not contain any load when visiting a satellite.

i∈v 1,2 X (i, s, k) = j∈v 1,2 X (s, j, k) ∀s ∈ sat , , k ∈ F 1,2 (54.7) W (s, k) ≤ j∈v 1,2 (X (s, j, k)) * WM AX ∀s ∈ sat , , k ∈ F 1,2 (54.8) U (i, k) ≤ M * (1 -X (i, e f , k)) ∀e f ∈ E f , i ∈ v 1,2 , k ∈ F 1,2
(54.9)

U (i, k) ≤ M * (1 -X (i, s, k)) ∀i ∈ v 2 , k ∈ F 2 , s ∈ sat , (54.10)
Then, C2P stations constraints were defined (i.e. opening stations, allocation of customers and service at the opened ones). Subsequently (54.11) limit to one station per client and (54.12) guarantee that for opened stations, customers are within the coverage radius.

S(i, j ) * Ds(i, j ) ≤ R( j ) ∀i ∈ cloc, j ∈ c2 p (54.12)

Results

See Table 54.2.

Our model was implemented in GUROBI 10.0 coded in Python 3.8.8 and tested under an Intel(R) Core (TM) i7-VPRO CPU 2.3 GHz 16 GB RAM is efficient on instances of up to 25 nodes. For all the C2P clients, the problem is first defined as an assignment problem and then a routing problem for the opened stations, which are accessible at all times (no time windows) and create more flexibility for routing. Similarly, grey zone customers add further flexibility to the model, as they can be part of any vehicle itinerary. However, home-delivery clients reduce model flexibility as their arrival times must be within their time windows. It is relevant to highlight that the definition of the maximum waiting time allowed plays a determining factor in the correct synchronization of the vehicles at the satellites and, consequently, in the route calculation. As well, it is noted that if the model does not allow vehicles to arrive before the lower limit of the time window (which generates waiting times at the customer nodes) the results tend to increase the use of vehicles, to create longer waiting times at the satellites and, in certain cases, due to the inflexibility of the problem, to end up in a non-feasible zone.

Conclusion and Future Research

Our model aims to offer a last mile delivery scheme, where the costs associated with waiting times at customer nodes and satellites are minimized through synchronization at the satellites and waiting costs reflected in the objective. Moreover, our model includes delivery options through the integration of C2P stations in addition to traditional home deliveries under time windows. Our mathematical formulation is efficient for small size instances, however due to the NP-hard nature of the problem, it is necessary to build approximate algorithms for solving large size and real-city instances to compare the performance of our model under conditions closer to those faced by real last mile delivery services.

  , j ) = 1 ∀i ∈ cloc (54.11)

Table 54

 54

	.1 Notation
	Sets
	v all	Set of all nodes
	E/E f	Set of the real and the cloned depots of each Echelon
	chd/chd 0,1,2	Set of all home delivery clients/set of home delivery clients in the Grey
		Zone (0), 1st Echelon (1) and 2nd Echelon (2)
	c2p	Set of all possible locations for lockers and satellites
	sat/sat ,	Set of real and cloned satellites each physical satellite (sat) is duplicated n
		times (once for each 2nd and Grey Zone customer).
	loc	Set of possible location for lockers
	cloc/cloc 1,2	Set of all c2p clients/set of c2p clients from the 1st (1) and 2nd Echelon (2)
	v 0,1,2	Set of nodes that can be visited by all vehicles (0), only by 1st Echelon
		vehicles (1) and only by 2nd Echelon vehicles (2)
	F/F 1;2	Fleet of all vehicles/fleet of 1st (1) and 2nd Echelon (2) vehicles
	Parameters
	M	Big number/sum of all customer demands
	TMAX/WMAX	Maximum route duration/"(TMAX/10)" as a value of the waiting time
		allowed at satellites

Time/Dis(i,j)

Travel time and distance from node (i) to node (j)

Table 54 .

 54

		F 1
	Computational results	chd 1 chd 2
	2
		chd 0