Mona Ibrahim
email: mona.ibrahim@univ-fcomte.fr

Nadia Yousfi Steiner
email: nadia.steiner@univ-fcomte.fr

Samir Jemei
email: samir.jemei@univ-fcomte.fr

Daniel Hissel
email: daniel.hissel@univ-fcomte.fr

Wavelets-based approach for online Fuel Cells Remaining Useful lifetime Prediction

Keywords: fuel cell prognostics, data-driven models, discrete wavelet transform, RUL prediction

Prognostics and Health Management (PHM) techniques for Proton Exchange Membrane Fuel Cell (PEMFC) systems are of great importance for increasing their reliability and sustainability. PEMFC systems suffer from relatively poor longterm performance and durability, and prediction and prognosis can give early indications about when components should be fixed or replaced. Prognostics modelling needs to take account of a number of phenomena, including degradation mechanisms that are not easily measured. A number of works are currently investigating PHM in fuel cell systems, as well as the problem of estimating Remaining Useful Lifetime (RUL). Any reduction in the volume of data required for making predictions is clearly advantageous. In this work, a univariate prognostic approach based on signal processing, namely Discrete Wavelet Transform (DWT) is proposed. The proposed approach aims at achieving an online prognostic for PEMFC systems. DWT is first introduced, and then the predictions are built using the power signals of two different PEMFC stacks in two different scenarios, namely static and dynamic operating conditions. Results show that the method is reliable for online prediction of power, with prediction errors less than 3%.

I. INTRODUCTION

n recent years PEMFCs have been seen as an environmentally friendly technology because of their clean chemical reactions (they use hydrogen to produce electricity, heat and water), and their efficiency. They can replace internal combustion engines in vehicles [START_REF] Ramirez-Murillo | Energy management of a fuel cell serial-parallel hybrid system[END_REF] [START_REF] Barbir | PEM Fuel Theory and Practice[END_REF]. Fuel cells (FCs) have the potential to provide heat and electricity not only for vehicles, but also for buildings, including houses.

Unfortunately, in comparison to other technologies, they have low reliability and a short lifespan [START_REF] Ibrahim | Energy storage system-characteristics and comparisons[END_REF] [START_REF] Beer | Fuel cell condition monitoring using optimized broad band impedance spectroscopy[END_REF].

The reliability of the PEMFC stack is improving over time, but like any device involving chemical, physical, and mechanical processes, the lifetime of PEMFCs may be shortened by a number of direct and indirect factors such as the deterioration of materials, imperfect cell/stack design and assembly, inadequate operating conditions, and impurities and contaminants in incoming gases. Although a loss of performance over time is unavoidable, the rate of degradation can be minimized through the use of effective diagnostic and prognostic tools that encapsulate an understanding of how and why degradation and failures occur. Diagnostic [START_REF] Li | Diagnosis for PEMFC systems: a data-driven approach with the capabilities of online adaptation and novel fault detection[END_REF]- [START_REF] Debanjak | Fuel cell condition monitoring system based on interconnected DC-DC converter and voltage monitoring[END_REF] and prognostic methodologies for assessing the state of health of a fuel cell and predicting its remaining lifetime need to be designed and implemented, to bring about an improvement in system performance.

Fuel cell system prognostics involve building a reliable model that is able to predict how the representative parameters for degradation and ageing of the FC will change, and that can estimate Remaining Useful Lifetime (RUL). Prognostics have several objectives: besides the aim of avoiding degradation and even catastrophic failures, there are the aims of extending FC lifetime and availability, optimizing service, minimizing risk and reducing costs. Approaches to prognosis can be separated into three broad groups: Model-based [START_REF] Yu | An integrated approach to prognosis of hybrid systems with unknown mode changes[END_REF]- [START_REF] Swindeman | A comparison of creep models for nickel base alloy for advanced energy systems[END_REF], data-driven [START_REF] Baraldi | Model-based and data driven prognostics under different available information[END_REF]- [START_REF] Liao | Discovering prognostic features using genetic programming in remaining useful life prediction[END_REF] and hybrid approach [START_REF] Jouin | Prognostics and health management of PEMFC-State of the art and remaining challenges[END_REF].

Since fuel cell systems are complex, non-linear, multi-physics (fluidic, thermal, mechanical, and electro-chemical) and multiscale (time, space) systems, with strong interaction between parameters belonging to different subsystems, it is difficult to develop complete, precise models of operation, degradation and ageing. For this reason, data-driven prognostics would appear easier to implement, since they do not require an a priori deep understanding of all the ageing mechanisms affecting the stack.

Javed et al. [START_REF] Javed | Data-driven prognostics of proton exchange membrane fuel cell stack with constraints based summation-Wavelet extreme learning machine[END_REF] have pointed out that data-driven approaches for estimating the RUL of a FC system can themselves be classified into two groups. The first group corresponds to univariate-degradation-based models and the second group corresponds to direct RUL prediction models.

Silva et al. [START_REF] Silva | Proton exchange membrane fuel cell degradation prediction based on adaptive neuro-fuzzy inference systems[END_REF] developed a data-driven approach that they called an Adaptive Neuro-Fuzzy Inference System (ANFIS), for PEMFC stack temperature prognostics, using data collected on the PEMFC stack.

In another work, Morando et al. [START_REF] Morando | Fuel cells prognostics using state network[END_REF] applied a data-driven approach, Echo State Network (ESN), for PEMFC stack prognostics.

In one recent paper [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF] the authors developed an approach based on a particle filtering for estimating the RUL of PEMFC under two operating conditions: steady state and dynamic.

Finally, Javed et al [START_REF] Javed | Enabling health monitoring approach based on vibration data for accurate prognostics[END_REF] presented an approach for PEMFC stack prognostics using a constraint-based Summation-Wavelet Extreme Learning Machine. The data-driven methods developed in the literature have two major problems. First, they all need a large amount of data for the learning phase (exogenous and past data). Second, they all need to be trained (learning from examples) to make predictions, which means that if new data become available (new examples), the models need to be rebuilt.

To help overcome these problems, in this paper we propose a new data-driven prognosis approach for PEMFC. The model is univariate: unlike multivariate approaches [START_REF] Javed | Enabling health monitoring approach based on vibration data for accurate prognostics[END_REF] that need exogenous data, our method does not take exogenous data into account, and predictions can be made using only past observations of the power signal (univariate time series), meaning that fewer data and acquisition devices are required.

The model uses the past values of this signal up to a defined time t and predicts new values up to a certain time t+h; thus, the process of degradation is predicted and an online prognosis can be made. The proposed approach is mainly based on Wavelet I Transform (WT), which has been widely used for fuel cell monitoring [START_REF] Beniouioua | PEMFC stack voltage singularitiy measument and fault classification[END_REF]- [START_REF] Javed | Improving accuracy of long-term prognostics of PEMFC stack to estimare remaining useful life[END_REF], as well as many other applications [START_REF] Bouzida | Fault diagnosis in industrial induction machines through discrete wavelet transform[END_REF]- [START_REF] Bialasiewicz | Wavelet-Based Approach to Evaluation of Signal Integrity[END_REF].

The advantage of WT, specifically of Discrete Wavelet Transform (DWT), is its ability to represent a given time series (or signal) with a small number of coefficients (the approximation coefficients), and consequently only a small quantity of data can be sufficient for performing a prediction. The proposed approach consists of considering the output signal of a PEMFC and decomposing a segment of 168 hours (one week) using DWT. The power signals considered in this article represent the data from two long-duration ageing tests of 1000 hours, and therefore they include the processes of degradation that will be considered in the present work.

After decomposition, a set of approximations coefficients (see section II) can be obtained and can represent the original signal with a smaller number of observations (168 hours represented by 21 hours). Predictions are made on this new set using mathematical regression models: Auto-Regressive Integrated Moving Average 'ARIMA' and polynomial regression. The prediction horizon is set to 21 hours. Once the prediction is made, inverse DWT is applied to create the approximation signal of 168 hours using the obtained prediction (see section II). These 168 hours represent the prediction of degradation, and will be compared with the actual signal in terms of the mean square error. Thereafter, this operation is repeated for the next 5 weeks in order to complete the online prognosis. The quantity of delivered power (with constant current) is an indicator of the RUL, since the performance of the considered PEMFC ceases to be satisfactory once the power drops by about 5.5% from the initial value (this value is given here as an example, and may be adapted according to the constraints of the application).

The novelty of the proposed approach is that unlike the approaches developed in the literature, it is a univariate model without a training part; that is, the model does not learn the relationship between the observations by training (as black box systems do, for example), but only by attributing a suitable analytical model. Moreover, a long-term prediction (168 hours) can be achieved by predicting only 21 hours using DWT and inverse DWT.

This paper is organized as follows: In section II the description of the wavelet transform and the model description are given and in section III the collection of experimental data is described. The analysis of results and discussion are detailed in section IV, and finally section V is a general conclusion of the present work.

II. PEMFC prognostics method: the mathematical framework

Discrete Wavelet Transform (DWT) is an effective method for non-stationary signal filtering and processing. DWT has been used in a variety of fields including medicine, energy management, image processing and computer vision. The greatest advantage of DWT is localization in both time and frequency [START_REF] Misiti | Wavelets and their applications[END_REF]. However, signals used for prognostics (voltage, power, etc…) are also non-stationary, and so DWT would appear to be a suitable method for analyzing these signals too.

A. Discrete wavelet transform

Discrete Wavelet Transform (DWT) is a method commonly used in signal processing. Given a time series (or a signal) Y (t), DWT allows an analysis of Y (t) in different scales over time. A wavelet φ is an oscillating, well-localized function having a finite vanishing moment [START_REF] Misiti | Wavelets and their applications[END_REF]- [START_REF] Zhu | Online condition monitoring in micro milling: a force waveform shape analysis appraoch[END_REF].

Using a function φ, with another wavelet ω, an orthogonal basis of L2 (the set of signals with finite energy) can be constructed (in this case, φ is termed the mother wavelet and ω the father wavelet). Using these two functions, the wavelet basis is defined by: 𝑊𝐵 = {𝜑 𝑗,𝑘 (𝑡) = 𝜑(2 -𝑗 𝑡 -𝑘), 𝜔 𝑗,𝑘 (𝑡) = 𝜔(2 -𝑗 𝑡 -𝑘) , 𝑗, 𝑘 ∈ ℤ} [START_REF] Ramirez-Murillo | Energy management of a fuel cell serial-parallel hybrid system[END_REF] where k is the translation parameter (in time), j is the scale (an integer indicating how the wavelet is stretched or compressed; see Fig. 1). The projection of the signal Y over the basis WB is called the discrete wavelet transform.

Choosing scale j, the DWT of Y can be written as follows:

𝑌(𝑡) = ∑ 2 -𝑗 2 𝑘∈ℤ 𝑐 𝑗 (𝑘)𝜑(2 -𝑗 𝑡 -𝑘) + ∑ ∑ 2 -𝑗 2 𝑘∈ℤ 𝑗 𝑖=1
𝑑 𝑖 (𝑘)𝜔(2 -𝑗 𝑡 -𝑘) (2) The first series c j (k) is known as the approximation coefficients and can represent the trend of Y, while d i (k) are known as the detail coefficients. Let C be the set of the approximation coefficients c j (k).

Equation (2) includes two operations: a down sampling (data divided by 2) and a convolution over the whole time interval [START_REF] Sheng | Wavelet Transform, the wavelet transforms and applications handbook[END_REF]. These characteristics are similar to a high pass and low pass filtering. Therefore, the wavelets (φ and ω) can be linked to low and high pass filters respectively [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF].

For each level j, the length of the original signal is divided by 2 as a result of the down sampling, and therefore the length of the approximation coefficients is equal to the length of the original signal divided by 2 j .

The DWT consists of two steps: analysis (or decompositions) where the detail and approximation coefficients are obtained by the filtering processes, and synthesis (or reconstruction) where the coefficient sets can be represented by signals having the same length of the original signal Y, and are obtained by a up sampling and a convolution [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF]. The process can be illustrated as shown in Fig. 2.

Inverse Discrete Wavelet Transform (IDWT) allows a reconstruction of the original signal using the set of approximation and detail coefficients, by an inverted operation of DWT.

It should be noted that the obtained signal A (signal constructed using the approximation coefficients) has the same shape as Y: in other words, it can describe the trend of Y, scale (j) scale (j) scale (j) without noises, and has the same length as Y.

translation (k) translation (k) translation (k)

B. Proposed PEMFC prognostic approach

The two fuel cell stacks considered in this study are from the same manufacturer and have the same initial characteristics. The first stack, FCd, was operated under dynamic current testing conditions, while the second stack, FCs, was operated in stationary regime under more or less nominal conditions. Both stacks were operated 24 hours a day over about 1000 hours.

The predictions were made using the power signals of both stacks FCd and FCs. The power was calculated using the formula:

𝑃(𝑡) = 𝑈(𝑡). 𝑈(𝑡) (3)
where U and I are respectively the given voltage and current of the stack. The voltage and the current are functions of time; therefore, the power P is also a function of time and can be considered as a signal or time series. P (t) is the observation at time t.

As we saw above, by applying the DWT of level j, the original signal of length L can be divided into j + 1 sets: the set of approximation coefficients of length L/2 j and j sets of detail coefficients.

It is important to note the difference between the approximation coefficients set C and the approximation signal (A). The signal A is constructed by using the set C and has a length of 2 j × l where l is the length of the set C, and this is obtained by the up sampling operation in the DWT (see Fig. 3).

The process that gives rise to the fuel cell RUL prediction can be summarized as follows:

The first step is to construct the prediction model, with the first week (168 hours) being taken as a learning window. The second step is to use this window to perform the prediction for the following week (the following 168 hours). Since in real time the data (or observations) are updated dynamically as the process unfolds, the prediction model uses the second available week (168 hours) to predict the third week (168 hours) and so on, up to 5 weeks.

In order to perform the RUL prediction, four prediction models were used. This allows us to compare the prediction results (first and the second models) and to show the impact of the wavelet transform on the two considered models (the third and the fourth models): 1) First model: A polynomial regression model uses the 168 power signal observations to predict the next 168 hours.

𝑦 𝑡 = 𝑐 + 𝛼 1 𝑦 𝑡-1 + 𝛼 2 𝑦 𝑡-2 + ⋯ + 𝛼 𝑝 𝑦 𝑡-𝑝 + 𝜃 1 𝑢 𝑡-1 + 𝜃 2 𝑢 𝑡-2 + ⋯ + 𝜃 𝑞 𝑢 𝑡-𝑞 + 𝑢 𝑡 (4)
where y t is a univariate time series, p is the autoregressive polynomial order, q is the moving average order, α 1 , α 2 , … , α p and θ 1 , θ 2 , … , θ q are the parameters (or coefficients) of the model, u t is a white noise.

It should be noted that d in ARIMA(p,d,q) is the differentiation order when the time series is not stationary.

3) Third model: A 3-level DWT using the 'db3' (Daubechie's 3) wavelet (these choices are detailed in section IV) is applied on the learning window; thus, the length l of the approximation is l=168/2 3 =21.

A new polynomial regression model is constructed using these 21 observations, and the model predicts only the next 21 points. Thus, a new set of approximation coefficients D is obtained. D is formed by the old set C plus the 21 prediction points. Thus the length of D is equal to 42.

The approximation signal A is now constructed using the set D. The length of A is 42 × 2 3 = 336, which shows that 168 hours are predicted using a prediction of only 21 points. Fig. 4 illustrates this procedure: it shows a set D that contains the set C of Fig. 3. We have the prediction of 21 points, and the reconstruction of the approximation signal A from the set D.

III. Experimental data

As part of the CNRS CO-CONPAC project [START_REF] Wadhame | Impact of power converter current ripple on the durability of a fuel cell stack[END_REF], durability tests are currently being carried out on two similar 600W 5-cell Their maximal current density is 1 A/cm². The test bench (Fig. 5) is adapted for PEM fuel cells with a power of up to 1 kW. A number of physical parameters involved in the stack can be controlled and measured in order to create the fuel cells' precise operating conditions as accurately as possible:

stack temperature, gas flows (from 0 to 100 l/min for air; and from 0 to 30 l/min for H2), fluid hygrometry rates (from 0 to 100%RH), air dew point temperature can be imposed; inlet and outlet flows, pressures (from 0 to 2 bars), temperatures (gas temperature from 20°C to 80°C), single cell voltages, stack voltage, current (from 0 to 300 A) can be monitored using a home-made interface developed with Labview™. Two independent boilers upstream of the stack provide gas humidification. Air and hydrogen go through their respective boilers before reaching the stack. The air boiler is heated in order to obtain the relative humidity desired, while the hydrogen boiler remains at room temperature. The temperature of the cell is controlled by a water-cooling system (with a cooling flow that can be set from 0 to 10L/min and a cooling temperature from 20°C to 80°C).

The current supplied by the battery is controlled by a TDI Dynaload-type active load. For the tests with high frequency ripple current, the load is driven by a low generator frequency. In normal operation, without current disturbance, the charge is conventionally driven from the management interface Labview.

A. Description of the tests carried out

Two long-term tests, each of 1000 hours, were carried out (Fig. 6 and Fig. 7).

1.

The ageing test was performed under dynamic condition

The FC stack durability test consists in operating an FC stack with identical characteristics under dynamic current loads and with the time target of 1000 hours. A ripple current of 70A with 7A oscillations at a frequency of 5 kHz was imposed on the FC in order to simulate the global operation effect of the high frequency power converter connected to the output of the FC stack. A complete characterization of the FC was performed every week (around 160 hours) in this order: polarization curve followed by electrochemical impedance spectroscopy (EIS) (the values are the same as for the second test). The characterizations were done at hours: 0, 35, 182, 343, 515, 666, 830, and 1016.

2.

The reference ageing test was performed on the second fuel cell stack (FCs). The durability experiment consists in operating the FC stack over 1000 hours in stationary conditions (a current of 70A is imposed). A complete characterization of the FC was performed every week. The duration of the tests was adapted to create the best fit with the first experiment. First, an electrochemical impedance spectroscopy (EIS) was performed only at 70A (0.70 A/cm²), in order to evaluate the state of the FC before it could be modified by the polarization curve. Second, a polarization curve was performed: the stack and cells voltages were measured under a current ramp from 0 A/cm² to 1 A/cm² of 1000s. In order to keep stoichiometric factors constant, the air and hydrogen flows were reduced to a current of 20A. An electrochemical impedance spectroscopy (EIS) was then performed, with measures made in the following order: 70 A (0.70 A/cm²), 45 A (0.45 A/cm²) and 20 A (0.20 A/cm²). Between each measure, a stabilization period of 15 minutes ensured the stability of the parameters. The characterizations were done at hours: 0, 48, 185, 348, 151, 658, 823, and 991.

B. Results / Ageing curves

Temporal curves of the voltage were analyzed in order to estimate the degradation rate. An illustration is given in Fig. 8, which depicts the voltage drop and the corresponding approximation using a least squares method. R² (a coefficient of determination between 0 and 1 that reflects how well the regression fits the data -the closer it is to 1, the better the fit) is higher for the second test, as a result of the absence of oscillations and discontinuities as the cells' voltage changes over time. A similar analysis (estimation of the degradation rate) was performed using polarization curves (see Fig. 9 for an illustration). In the two experiments, we can see that the fifth cell presents a higher degradation rate; this is due to its position in the stack.

The experiment was run for around 1000 hours under stationary conditions, and the average stack voltage decay rate was shown to be 153 μV/h per stack and 29.6 μV/h per cell. However, the voltage decay rate observed in dynamic conditions is slightly higher (less than 1% of average voltage) at 172 μV/h per stack and 33.4 μV/h per cell.

In order to implement the proposed approach, a control system was installed based on an Autobox™ (Dspace) using MATLAB/SIMULINK, enabling a direct connection between the model and the system.

IV. Results and discussion

The prediction results corresponding to the different proposed models are presented in Fig. 10 and Fig. 11.

Tables (I) and (II) show how the degradations (estimated using the mean values) progressed over time. These degradations are measured with respect to the initial values of the considered power signals.

The online degradation predictions are made as follows: The first week (168 hours) is used to construct and learn the models. The predictions are performed over the following week (the second week). In other words, the prediction horizon is 168 hours.

The models are updated using the new available data for the second week, and the predictions are made over the third week; and so on, up to five weeks.

The model simulations and predictions were done using MATLAB software A list of existing programmed MATLAB functions can be found in the econometric and wavelet toolboxes, to be used for time series modelling, prediction and decompositions.

Results are analyzed in the following two paragraphs. The models are compared to show how DWT was able to improve the predictions. In the first paragraph, models 1 and 3 (without DWT) are compared with models 2 and 4, where the DWT is used (see section II), for the first and second stacks FCs and FCd. In the second paragraph, the results for FCs and FCd are compared.

A. Analysis of the results: comparison of the models

In order to compare the prediction results, the first and second models are compared with the third and fourth models respectively.

Regarding the first model (regression without DWT), a polynomial of degree 1 was attributed to the learning window: this polynomial fits the data better than higher polynomial degrees. For the DWT part in the third model, a 3-level decomposition and construction and the Daubechie's 3 wavelet (db3) were used for both the decomposition and the reconstruction operations. These choices are justified as follows:

Since the division of 168 hours (learning part) by 2 3 gives 21 hours, this is the nearest obtained value to 24 hours (one day), so a general trend of the power during one day can be observed. The mother wavelet is db3. Sym3 would give the same results, since db3 and sym3 have exactly the same filter coefficients [START_REF] Misiti | Wavelets and their applications[END_REF]. Table (III) gives values of the mean square error defined by equation ( 5) between the approximation and the original signal, for the same level of decomposition using different mother wavelets: Morlet, Coif (Coiflet), Meyer, Sym (Symlet), and Db (Daubechie's). The numbers in the wavelets Coif2, 3, 4 and 5, Sym 2, 3, 4 and 5 and Db2, 3, 4, and 5 correspond to the numbers of terms of the corresponding wavelet [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF]. The values in table III represent the errors (obtained using equation ( 5)) between the signal of approximation (reconstructed using the approximation coefficients) and the original signal. These values are provided in order to help finding the best wavelet to perform the algorithm. The wavelet that gives the smallest error, is the then chosen.

𝐸 = 1 𝑁 ∑ (𝑃(𝑡) -𝑎 𝑡 ) 2 𝑁 𝑖=1 (5)
where N is the length of the approximation signal of level 3, 𝑎 𝑡 is the t th approximation signal and 𝑃(𝑡) is the actual power. The chosen wavelet is the one gives the smallest error for the approximation signal, which is db3.

In Fig. 10(a) and 10(c), it can be seen that for the third model, the predictions are closer to the actual observations. Although they do not have exactly the same variations of the original signals, their behavior can nevertheless correspond: for much of the time they can approximate the mean value of the signals, and occasionally be very similar to them (see table I), which can provide solid indications about processes of degradation. Table (IV) shows the prediction errors.

The error used in this study is the relative error having the formulations:

𝑅𝑀𝑆𝐸(𝑒) = √ ∑ (𝑃(𝑡)-𝑃 ̂(𝑡)) 2
𝑛 𝑡=1 𝑛 [START_REF] Kunusch | modeling , diagnosis and control of fuel-cell-based technologies and theur integrations in smart-grids and automotive systems[END_REF] and max(𝑒) = 𝑚𝑎𝑥 | 𝑃(𝑡)-𝑃 ̂(𝑡) 𝑃(𝑡) | [START_REF] Debanjak | Fuel cell condition monitoring system based on interconnected DC-DC converter and voltage monitoring[END_REF] where 𝑃 ̂(𝑡) is the prediction of P at time t.

Tables IV to VI show that the predictions made using the third model have the fewest errors. Regarding the comparison between the second and the fourth models:

The ARIMA models that best fit the observations are ARIMA(5,1,0), ARIMA(5,0,1) and ARIMA(5,0,0). The model updates its coefficients when new actual values of a new week are made available. As shown in Fig. 10(b) and 10(d), DWT can improve the prediction for the ARIMA models. Table (I) shows that the predictions can also approximate the mean of the actual observations, and in table (IV) the comparison of errors proves again that DWT improves the quality of prediction. The comparisons of models for the stack FCd (dynamic conditions) can be seen in Fig. 11, and in tables (II) and (VI).

For Fig. 11 (a), the polynomial regression fails to capture the general behaviors of the degradations obtained by the first model. However, Fig. 11(c) shows the improvement obtained when DWT is added to the polynomial regression (third model); here the polynomial is of degree 2. It can clearly be seen that with DWT the model captures the process of degradation. Tables (I) and (IV) give numerically the impact of DWT in reducing the prediction errors.

Here too, the impact of DWT is demonstrated. The situation for Fig. 11 Finally, as shown in tables (I) to (VI),the errors obtained by the ARIMA-DWT model are less than the regression-DWT, showing thus that this is the best model (among the four models considered in this paper) to implement a prognostic strategy for a fuel cell.

B. Analysis of the results: comparison of the stacks

The main difference between stacks FCs and FCd is that stack FCd is subjected to dynamic operating conditions (see section II) that affect the prediction results.

Let us now look at Fig. 10(d) and 11(d). In the case of stack FCs, the model shows more stability in approximating the observations. But this is to be expected, since the trend of the voltage signal of FCd (see Fig. 8) drops linearly and no ripple variations are added to this stack; while the trend of the voltage for FCs is variable, since dynamic ripples are added to the stack. Nevertheless, the predictions can be considered as acceptable for both stacks, according to the small percentage errors given in table (VI).

C. RUL Estimation:

Since the best prediction results are obtained using the fourth model, the RUL is estimated according to this model. A flowchart of the algorithm represented by the fourth model can be found in Fig. 12. In this work the RUL is defined as the time (in hours) until power drops by 5.5% in relation to its initial value. As for the stack FCd, the results are presented in Fig. 13(b). The initial value of the power is also about 236 W and after a drop of 5.5% it is therefore about 223 W, indicating an actual RUL of about 400 hours (the third week). The RUL estimated using the fourth model was found to be 469 hours (also in the third week).

We conclude that model 4 can give a quite robust estimation of the RUL, close to the real figure, for both stacks under different operating conditions, static and dynamic.

D. Discussions

The obtained results are very promising. A number of strong points of the proposed method can be summarized as follows:

• Its rapid execution time for both the learning and the prediction parts (an estimation of the order of seconds) makes it superior to other data-driven algorithm, that is to say artificial neural network algorithms that have a high time cost. • It has no training part and thus does not include any black box systems that are often difficult to adapt to any change in data. Instead, it uses analytical equations. • The model requires fewer data (only 168 observations) for learning, and no exogenous data are needed to perform the prediction for both stacks. A reduced reliance on measurement devices, means reduced instrumentation cost and reduced testing time.

• It provides a close approximation of the actual observations, as shown in table (VI) (no more than 2.81% error). Moreover, it is clear that, even under dynamic operations, the predictions include variations that are similar to the actual variations (Fig. 11(d)). • The model is able to estimate RUL, since the power decrease has direct impact on the RUL. In the case of the stacks considered in this work, the RUL is estimated as the time until the power drops by 5.5%. The proposed model could nevertheless be improved by adding some degradation factors such as the use of materials and the occurrence of drying and flooding faults in the fuel cell: these two faults occur successively inside the fuel cell and have negative impact on the performance of the cell [START_REF] Gerard | Ripple current effects on pemfc aging test by experimental and modeling[END_REF], leading to drops in the power signals, indicating a loss of fuel cell efficiency.

V. Conclusion

An approach based on DWT for estimating the RUL for a stack of PEMFCs is presented in this paper. The novelty of the presented techniques is the use of DWT for estimating the RUL, and the use of univariate models without training. The main advantage, in comparison to similar works on this topic, is the time series domain analysis that uses past information from a univariate time series (the power signal). Its reliance on measurement devices and on data are therefore reduced. This approach can be used for general cases, since the mathematical models do not take into account the variations that are to be found in exogenous data, under different operating conditions. It can be coded easily in any computer language, and thus easily implemented in real life applications.

The proposed method was simulated in MATLAB/SIMULINK, and compared with other models;

the results showed that with DWT, the prediction is more accurate and faster. It may be interesting to compare the DWT with other filters, such as unscented particle filters and adaptive neuro-fuzzy inference systems that have also been used for estimating RUL in electrochemical batteries.

Fig. 1 .

 1 Fig. 1. Translation and scaling of the wavelet transform

Fig. 2 .

 2 Fig. 2. The approximation coefficients and approximation signal obtained by DWT for a given time series.

Fig. 3 .

 3 Fig. 3. The approximation coefficients and approximation signal obtained byDWT for a given time series.

2 )

 2 Second model: An Auto-Regressive Integrated Moving Average (ARIMA) model functions in the same way as the first model. The ARIMA model, or, more precisely, ARIMA (p,d,q) [41][42] is formulated as follows: Approximation signal : the signal A Reconstruction Decomposition

4 )

 4 Fourth model: An ARIMA model is used to model the set C used in the third model; the set D is obtained by adding the prediction of 21 points to the set C. The prediction also uses an ARIMA model, and is compared with the prediction by the polynomial regression of the third model.

Fig. 4 .

 4 Fig. 4. The approximation coefficients, the set D of predicted coefficients and approximation signal obtained by DWT for a given time series. stacks (PEMFCs are assembled at the FCLAB and comprise commercial membranes, diffusion layers and machined flow distribution plates), each cell having an active surface of 100 cm². The first stack FCd is operated under dynamic current operating conditions (ripple current at 5 kHz superimposed on a constant currentpeak-to-peak value of 10% of constant current). A second fuel cell stack FCs is operated in stationary regime at more or less nominal conditions over 1000 hours; this second test acts as a reference for the first test. A series of observations are carried out once a week (around every 160 hours) according to the same protocol each time: a polarization curve test, and global historic curves and Electrochemical Impedance Spectroscopy (EIS) measurements. Relationships between ripple current and fuel cell performance, such as power loss and degradation, are investigated. The nominal current density of the fuel cells is 0.70 A/cm².Their maximal current density is 1 A/cm². The test bench (Fig.5) is adapted for PEM fuel cells with a power of up to 1 kW.

Fig. 5 :

 5 Fig. 5: View of the 1 kW FC test bench.

Fig. 6 :

 6 Fig. 6: 1st long-term test (FCd): 1000 h with high-frequency current oscillations

Fig. 8 :

 8 Fig. 8: Comparison of cell average voltage (Red : FCd with current oscillations ; Blue : FCs without current oscillations)

Fig. 9 :

 9 Fig. 9: Comparison of cell average voltage (polarization curves at I=70 A)

Fig. 11

 11 Fig. 11 Prediction results of the power signal for FCd (dynamic conditions) applying the first, second third and fourth models in (a), (b), (c) and (d) respectively. The red and the blue lines represent respectively the prediction and the actual signal

  (b) (second model) and Fig. 11(d) (fourth model) is similar: DWT is seen to improve the prediction with the ARIMA model.

Fig. 12 A

 12 Fig. 12 A flowchart of the fourth model.

Fig. 13 (

 13 Fig.13(a) shows the results for stack FCs. As the initial value of the power is about 236 W, the actual RUL is about 811 hours, corresponding to a drop of 5.5% to ~223 W (end of the fifth week). Using model 4, the predicted (or estimated) RUL for the same stack is 893 hours (beginning of the sixth week), which represents a difference of ~82 hours between the estimation and the actual RUL.As for the stack FCd, the results are presented in Fig.13(b). The initial value of the power is also about 236 W and after a drop of 5.5% it is therefore about 223 W, indicating an actual RUL of about 400 hours (the third week). The RUL estimated using the fourth model was found to be 469 hours (also in the third week).

Fig. 13 (

 13 Fig. 13 (a) RUL prediction of the stack FCs using the fourth model, and (b) RUL prediction of the stack FCd using the fourth model. The red line correspond to the prediction and the blue line corresponds to the actual signal


Jean-Baptiste Horel
email: jean-baptiste.horel@inria.fr

Christian Laugier
email: christian.laugier@inria.fr

Lina Marsso
email: lina.marsso@utoronto.ca

Radu Mateescu
email: radu.mateescu@inria.fr

Lucie Muller
email: lucie.muller@inria.fr

Anshul Paigwar
email: anshul.paigwar@inria.fr

Alessandro Renzaglia
email: alessandro.renzaglia@inria.fr

Wendelin Serwe
email: wendelin.serwe@inria.fr

Using Formal Conformance Testing to Generate Scenarios for Autonomous Vehicles

Keywords: Behavior trees, CARLA simulator, Formal methods, Input-output conformance, Scenario generation, Test purpose

Simulation, a common practice to evaluate autonomous vehicles, requires to specify realistic scenarios, in particular critical ones, occurring rarely and potentially dangerous to reproduce on the road. Such scenarios may be either generated randomly, or specified manually. Randomly generating scenarios is easy, but their relevance might be difficult to assess. Manually specified scenarios can focus on a given feature, but their design might be difficult and time-consuming, especially to achieve satisfactory coverage. In this work, we propose an automatic approach to generate a large number of relevant critical scenarios for autonomous driving simulators. The approach is based on the generation of behavioral conformance tests from a formal model (specifying the ground truth configuration with the range of vehicle behaviors) and a test purpose (specifying the critical feature to focus on). The obtained abstract test cases cover, by construction, all possible executions exercising a given feature, and can be automatically translated into the inputs of autonomous driving simulators. We illustrate our approach by generating thousands of behavior trees for the CARLA simulator for several realistic configurations.

I. INTRODUCTION

AV (Autonomous Vehicles) are complex and safety critical systems. Although their main objective is to propose a safer alternative to human-driven cars, in the last years autonomous cars were involved in several accidents, some of them even fatal (e.g., the ones of a Uber car in 2018 and a Tesla in A part of the work has been performed in the project ArchitectECA2030 that has been accepted for funding within the Electronic Components and Systems for European Leadership Joint Undertaking in collaboration with the European Union's H2020 Framework Programme (H2020/2014-2020) and National Authorities, under grant agreement No. 877539. A part of this work has been supported by the PRISSMA project, co-financed by the French Grand Défi on Trustworthy AI for Industry. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

* Institute of Engineering Univ. Grenoble Alpes 2019 [START_REF] Boudette | It Happened So Fast[END_REF] 1 ). To prevent these accidents, all the components of an autonomous vehicle-which include both the sensors (e.g. LiDAR, cameras, GPS, etc.) and the algorithms that are behind the perception, decision making and control modules-need to be thoroughly tested to ensure they handle critical situations even better than human drivers. To avoid taking unnecessary risks, and because such critical situations are unlikely to happen-at least frequently enough-in real road traffic situations, a common and essential [START_REF]for the Harmonization of Vehicle Regulations[END_REF] practice is to reproduce these critical situations in a simulator, such as CARLA [START_REF] Dosovitskiy | CARLA: An open urban driving simulator[END_REF]. However, these scenarios to drive a simulator need to be specified. Two techniques are used currently: constrained random generation and manual specification [START_REF] Riedmaier | Survey on scenario-based safety assessment of automated vehicles[END_REF]. Randomly generated scenarios can be easily produced, but their relevance might be difficult to assess, since they can present a high level of redundancy, which is hard to detect and strongly limits their coverage [START_REF] Riedmaier | Survey on scenario-based safety assessment of automated vehicles[END_REF]. On the other hand, manually specifying a large number of scenarios is extremely time consuming and a satisfactory coverage of all possible situations is hardly achievable.

In this work we propose to apply formal methods to automatically generate scenarios, which are guaranteed to be relevant for testing AV's behaviour in a particular situation (e.g., collision, near miss, etc). More precisely, we will use a conformance testing tool to generate scenarios from a formal model and a test purpose characterizing the situation. The formal model is specialized to a given configuration, which includes a scene map and several actors with their initial positions and constraints on their trajectories. Precise trajectories of the actors will be automatically induced by the generated scenarios. In general, each test purpose will yield several scenarios, with guarantees to cover all relevant variations of the behavior related to the test purpose. These scenarios are then automatically transformed to be used as input for a driving simulator. To ensure the generation of relevant and critical scenarios, test purposes (e.g., reaching a collision) and test configurations can be defined based on critical situations emerging from road accident data [START_REF] Ding | Learning to collide: An adaptive safety-critical scenarios generating method[END_REF]. We illustrate our approach with CARLA by providing a method to translate the scenarios into behavior trees. Our approach is evaluated on ten configurations, involving three scene maps (T-crossing, highway, and X-crossing) and various actors, for which we generated several scenarios featuring collisions of the AV with other actors, near-misses of such collisions, and arrivals at the destination.

The rest of the paper is organized as follows. Section II compares our approach to existing ones for scenario generation in driving simulators. Section III gives a bird's eye view of our approach. Section IV illustrates the application of our approach on several configurations. Finally, Section V gives some concluding remarks and future work directions.

II. RELATED WORK

Although AV scenario generation approaches have been intensively studied [START_REF] Riedmaier | Survey on scenario-based safety assessment of automated vehicles[END_REF], current methods never test the complete ODD (Operational Design Domain) of the AV, but rather restrict scenarios to specific configurations, such as highway overtaking or an urban intersection. Hence, the testing methods can neither verify the full safety of the AVs nor be used at industrial scale, as they do not guarantee a complete coverage of the ODD.

A common way to create interesting and critical AV scenarios is to start from an existing abstract and parameterized scenario, from which several AV scenarios can be generated by choosing the values of the parameters. As not all of them are critical and interesting to run, several studies adopt optimization methods to find critical sets of parameters. The three following approaches use stochastic optimization methods to search the parameters values leading to the corner cases of an abstract scenario. The Matlab toolbox S-TaLiRo implements several optimization methods used in [START_REF] Tuncali | Utilizing s-taliro as an automatic test generation framework for autonomous vehicles[END_REF] to generate scenarios causing the motion controller of an AV to fail in an overtaking abstract scenario and generating a collision. Bayesian optimization is used by [START_REF] Gangopadhyay | Identification of test cases for automated driving systems using bayesian optimization[END_REF] to generate AV scenarios from a parameterized abstract scenario. In particular, to find the values of these parameters that cause collision between a car and a pedestrian, the optimization method also exploits the feedback from the simulated execution of the previously generated scenario. An AEB (Automatic Emergency Braking) system is tested in [START_REF] Khastgir | Test scenario generation for driving simulators using constrained randomization technique[END_REF] using constrained randomization techniques to generate AV scenarios from an abstract scenario and then converge to critical AV scenarios by reducing the constraints. Our approach generates the actors behaviors (like car trajectories), thus we do not use an abstract scenario with predefined and parameterized behaviors as input. Additionally, our approach can generate several different scenarios covering better the many possible outcomes of an abstract scenario.

[1] and [START_REF] Klischat | Generating critical test scenarios for automated vehicles with evolutionary algorithms[END_REF] relate the criticality of a scenario to the dimension of the solution space, defined as the control space that does not lead to a collision. The solution space is computed using reachability analysis. [START_REF] Klischat | Generating critical test scenarios for automated vehicles with evolutionary algorithms[END_REF] uses an evolutionary algorithm to find the best parameters (speeds, initial positions, etc.) of an existing scenario to minimize the solution space, while [START_REF] Althoff | Automatic generation of safety-critical test scenarios for collision avoidance of road vehicles[END_REF] uses optimization techniques by modeling the generation problem as a quadratic problem with constraints. In contrast, our approach does not start with an existing parameterized scenario to then increase its criticality. We can generate this scenario based on an abstract scenario, and one could reuse our generated scenarios as an input for these approaches. [START_REF] Li | Ontology-based test generation for automated and autonomous driving functions[END_REF] proposes a different method to test an AEB system by defining an ontology based on the ODD of the AEB system, and then by using a combinatorial test generation approach to create AV scenarios from the ontology. A machine learning-based approach is also used to increase the criticality of the scenarios by exploiting the result of the already simulated scenarios. In contrast, our approach is used on several abstract scenarios while the ontology is associated to one abstract scenario of AEB testing. In addition, the complexity of the scenarios depends on the number of parameters and their value range in the ontology, and so this approach depends more on heavily parameterized scenarios than our approach.

Machine learning is also used to create realistic and critical scenarios. In [START_REF] Ding | Learning to collide: An adaptive safety-critical scenarios generating method[END_REF], a neural network generates safety-critical AV scenarios from an abstract scenario of an urban intersection. The obtained scenarios are modeled as series of probability distributions, from which a final one is sampled and run on CARLA. [START_REF] Krajewski | Data-driven maneuver modeling using generative adversarial networks and variational autoencoders for safety validation of highly automated vehicles[END_REF] uses a generative adversarial network, trained on a dataset of vehicle trajectories on a highway, to generate realistic lane change trajectories. While learning-based generators are limited to what they learned, our approach does not have standard machine learning drawbacks, does not require long training, does not over fit the training data, and hence can generalize to new domains, i.e., we can include new configurations at inference time. Most importantly, our approach comes with guarantees that the generated scenarios are diverse (i.e, different trajectories), and all relevant to the AV's behavior in a given event (specified in the test purpose).

III. METHODOLOGY

Figure 1 gives an overview of the proposed flow. Its first input is a configuration defining the scene with its objects and their behavior, from which a formal model and a corresponding CARLA configuration are derived. The second input is a test purpose, describing the intent all test cases should focus on. From these inputs, we automatically compute a comprehensive test suite, and translate each generated test case into a behavior tree to drive the CARLA simulator. The subsequent sections present these steps, which are all fully automatic, besides the two steps represented by dotted lines, i.e., deriving a formal model and a CARLA configuration. To ease the latter two steps, we provide libraries factoring common parts.

A. Formal model of a scene

We focus on scenes with an ego vehicle, called car, moving around in the scene towards a goal or destination position, trying to avoid any collision with the obstacles. The car features sensors (e.g., camera, LiDAR, etc.) enabling it to perceive its environment. Obstacles represent the different hazards that may disturb the progression of the car. There are two kinds of obstacles: the static, fixed ones (e.g., buildings, trees, parked cars, etc.), and the dynamic, moving ones (e.g., cars, cyclists, pedestrians, etc.). The trajectories of dynamic obstacles can be predefined or (partially) left random to leave more room for varying behavior. Both kinds of obstacles can be transparent, in which case the sensors of the car can perceive the status of the cells behind the obstacle, so as to enable taking into account the difference between, e.g., a ball and a truck. Considering a discrete, two-dimensional projection, we represent the map of a scene as a grid composed of cells, which can have one of the three different values: free, occupied by the car, or occupied by an obstacle. This discrete modeling helps to keep the size of the model tractable (note that the grid resolution can be made finer at the price of increased computing resources). The grid represents the ground truth under ideal conditions and is updated at each move of one of the actors (the car or a moving obstacle). Figure 2a shows the grid of a road-crossing scene. The buildings delimiting the two roads are represented as immobile obstacles (cells marked "I") and the car is denoted by the cell marked "C". The scene also comprises two moving obstacles (cells marked "M"): another vehicle in front of the car, and a pedestrian crossing the trajectory of the two vehicles. The cell of the pedestrian is depicted in a lighter color than the one of the other vehicle to indicate that the pedestrian is considered transparent and does not obstruct the view on the cells behind. A screenshot of a corresponding CARLA simulation is shown in Figure 2b.

We consider a car equipped with a LiDAR. For each configuration of the ground truth, the model can compute the perception grid, corresponding to the ideal expected output of the perception algorithm fed with the car's LiDAR input. We model the perception as a grid centered around the car, indicating whether a cell is free, occupied, or unknown (hidden from view). For each occupied cell, the perception also contains information whether the cell was occupied before, enabling to detect car or obstacle moves.

To formally model the behavior of such a scene, we use LNT [START_REF] Champelovier | Reference Manual of the LNT to LOTOS Translator (Version 7.0)[END_REF], [START_REF] Garavel | From LOTOS to LNT[END_REF], a language equipped with a formal semantics rooted in concurrency theory and a syntax close to classical programming languages. Actors (car and moving obstacles) are represented as concurrent LNT processes that interact by means of multiway rendezvous [START_REF] Charlesworth | The Multiway Rendezvous[END_REF], [START_REF] Garavel | The Unheralded Value of the Multiway Rendezvous: Illustration with the Production Cell Benchmark[END_REF], enabling synchronization and communication. Each actor process has a local copy of the map, to ensure that its moves obey the laws of physics. The current status of the map is handled by a map manager process, in charge of maintaining a consistent ground truth. When moving, each actor sends its new position to the map manager, which updates the map and broadcasts it to all actors. This allows to represent moves to any adjacent cells (including diagonals) or even farther cells. The scenario terminates when either the car arrives at its destination, or a collision occurs: in these cases, the model will perform one of the special actions ARRIVAL or COLLISION and then stop immediately. This implies that an execution sequence of the model will never contain both, a collision and the arrival of the car at its destination.

LNT has an interleaving semantics, in which actions are considered atomic and two different actions cannot be observed simultaneously. To connect to a simulator where moves take some time and can thus happen concurrently, we introduce a notion of discrete time and add a scheduler process, which generates special TICK actions and ensures that all actors move at most once between two ticks. Thus, all moves between two TICK actions can be executed in parallel in the simulator. The scheduler also enforces a fixed order of execution of the actors between two TICK actions, which prunes redundant interleaving of actor executions and makes the subsequent analyses more efficient.

To facilitate the handling of various configurations, the LNT model of a configuration is split in two different parts: a generic part, defining types, functions, and processes common to all configurations, and a particular part, defining the constants characterizing the considered configuration. The latter part defines in particular the two dimensions of the rectangular grid, the initial grid (indicating any static obstacles), the number of obstacles, as well as the initial position and behavior of all actors (car and obstacles).

The behavior of the car is defined as a sequence of actions, specified by the speed and direction of move. When all actions of the sequence have been performed, the car has arrived at its destination. The behavior of an obstacle is defined in the same way, with three extensions. First, an obstacle may choose not to move (between two TICK actions). Second, an obstacle may perform a random move, where the direction is chosen arbitrarily (respecting the physical constraints of the scene and avoiding collisions). Third, rather than stopping, an obstacle may restart the execution of its behavior sequence. In the current version of the model, the actors occupy a single cell and can turn as they wish. In principle, it is possible to refine the model and add constraints to represent rules of the dynamics (e.g., turning in at least x cells). The LNT model is compiled into an automaton (statetransition graph) using the compilers of the CADP toolbox [START_REF] Garavel | CADP 2011: A Toolbox for the Construction and Analysis of Distributed Processes[END_REF]. To generate tests from this automaton, we must specify which of its actions (i.e., transition labels) are considered as inputs or outputs. In our case, only the moves of the actors are inputs (to control the simulation), all the other actions being outputs (to observe the progress of the simulation), e.g., ground truth map, perception, detection of a collision, and arrival of the car. In the sequel, we will denote by "model" either the LNT description or its underlying automaton, since both represent the same behavior, albeit at different abstraction levels.

The model was validated using the CADP tools to ensure that it represents the desired scene and behaviors of the actors. Besides interactive simulation of the model (exploring back and forth the transition sequences from the initial state), we also checked the temporal logic property stating that a terminal state (either a collision, or the arrival of the car at its destination) will be eventually reached.

B. Computation of a Complete Test Graph

To obtain test scenarios focused on specific situations, we apply the conformance testing tool, which generates interaction scenarios with the simulator-in their simplest form, a sequence of transitions. As shown in Figure 1, we first extract a CTG (complete test graph) from the model and a TP (test purpose) using the TESTOR tool [START_REF] Marsso | TESTOR: A Modular Tool for On-the-Fly Conformance Test Case Generation[END_REF]. A TP is an automaton with special "ACCEPT" labels characterizing the states to be reached by the scenario, and a CTG is an automaton that contains all transition sequences leading to these states.

As an example, Figure 3 shows a TP to generate scenarios leading to a collision of the car with a pedestrian. This TP has two states and two transitions, requesting to reach (after an arbitrary number of transitions) a collision, represented by a transition labeled with "COLLISION !PEDESTRIAN".

When computing a CTG, only the transitions corresponding to a controllable input or observable output of the SUT (system under test, in our case the CARLA simulator) are necessary. Thus, we can hide-and reduce the model-all other transitions (e.g., the broadcast of the ground truth map) that are useful for validation, but irrelevant for test generation.

C. Extraction of test cases

In general, a CTG contains states for which several inputs can lead to a successful run. Thus, we apply the techniques of [START_REF] Marsso | Automated Transition Coverage in Behavioural Conformance Testing[END_REF] to extract a test suite, i.e., a set of TCs (test cases) covering all transitions of the CTG. Each TC is an automaton interacting with the SUT to drive it towards the accepting states specified by the TP. Thus, for a given model, using even a simple TP as the one shown in Figure 3, several different TCs (and hence, scenarios) can be generated. In our setting, all generated TCs are sequences, since from a given state there cannot be two different outputs.

D. Translation of test cases into behavior trees

The TCs extracted from the model and a TP are represented in an abstract form as automata, and must be transformed into a more concrete form to be used as simulation scenarios. This last step is dependent on the simulator considered; we describe below the procedure we developed for the CARLA simulator. A simulation scenario in CARLA begins by initializing the configuration (scene map and placement of actors) symmetrically w.r.t. the formal model, and then specifying the behaviors of the actors. To run a simulation scenario, we use (an extended version of) the scenario runner2 feature of CARLA, which takes as input a behavior tree and uses it to control the actor moves during the scenario.

Our procedure to translate a TC into a behavior tree consists of three steps. Firstly, the TC is converted to a textual form, more suitable for parsing. Given that in our setting a TC is always a sequence of transitions, its corresponding text file contains one action per line, ending with a PASS action that denotes the end of the sequence (i.e., the accepting state of the TP). Secondly, the text file is analyzed line-by-line using a shell script, which extracts the necessary information about the moves and new positions of the actors, and produces a JSON file containing an initialization part (declaration of the actors with their initial placement and orientation) and a behavior part describing the sequence of moves for each actor. Thirdly, this JSON file is fed as input to a Python program that parses it, builds a behavior tree in memory using the scenario runner primitives, and traverses this behavior tree to drive the simulation in CARLA by sending the commands associated to the tree nodes. Figure 4 shows an example of TC and Figure 5 shows its corresponding behavior tree. To make the simulation smooth and realistic, the moves of the actors between two consecutive TICK actions in the TC are executed simultaneously by encapsulating them in parallel nodes of the behavior tree. These parallel nodes are then encapsulated into one sequential node (the root of the behavior tree) representing the sequence of all simultaneous moves in the scenario.

IV. EXPERIMENTS

We illustrate our approach on several configurations, with different maps and obstacle behaviors. The experiments were performed on the troll cluster of the Grid'5000 platform, i.e., servers with an Intel Xeon Gold 5218 processor and 384 GB of RAM. The CARLA simulations were run on a laptop with an Intel i5-7440HQ, 16 GB of RAM, and an Nvidia Quadro P4000. The simulation was synchronized with CARLA's scenario runner to ensure that, for each simulation step, the simulation time was constant and the scenario runner completed its control loop. Fig. 4: Excerpt of a test case generated for model "X-crossing1" (see Fig. 2) and test purpose "collision pedestrian" (see Fig. 3) I gives statistics about the formal models (LNT descriptions and their underlying automata) of the considered configurations. For each model, the table indicates the number of obstacles, the number of LNT lines specific to the configuration (the generic part, common to all models, has 1130 LNT lines) and the size (number of states and transitions) of the minimized automaton. The execution time in seconds (resp., memory requirements in MB) corresponds to the sum (resp., peak) of both generation and minimization of the automaton.

All models except "highway" and "T-cross" use the map shown in Figure 2a. Models "X-cross-N " (N ∈ {1, 2, 3, 4}) focus on a collision between the ego car and a pedestrian, featuring also another moving car; differing in the initial positions and behavior, these configurations were designed as simple tests of the approach. In model "near miss" the itinerary of the car is crossed by different obstacles, leading to near misses (i.e., an obstacle passing very close to the car, but without colliding with it). Model "immobile car" features the car stopped in the middle of the crossroad, with various obstacles moving around, but without causing any collision. Model "X-cross-rand1" introduces additional variability by enabling random moves of the pedestrian and the other car; their precise moves in a TC will be chosen according to the TP. Model "X-cross-rand2" is similar to "X-cross-rand1", but features a different and longer trajectory of the ego car (11 moves instead of 7). Model "T-cross" uses a map with a T-shaped crossroad, with a collision between the ego car and another car ignoring priority rules. Model "highway" uses a rectangular map representing a long horizontal highway with the ego car on the middle lane, a faster car trying to overtake it, and a slower car ahead of it, but leaving the lane to let the ego car pass.

The difference in the model sizes can be explained by the varying number of obstacles and the number of moves of the actors. Models "immobile car" and "near miss" feature three, instead of two, obstacles, with a large number of moves (in total, 12 moves for "near miss" and 23 for "immobile car"). In particular, a random obstacle move amounts to possibly five transitions, corresponding to each direction it can go (up, down, right, left, or none). Also, as the model stops when the ego car has finished, the longer the car trajectory, the larger the size of the underlying automaton, as can be seen by comparing models "X-cross-rand1" and "X-cross-rand2".

We consider different test purposes for these models. Purpose "arrival" requires the car to complete its behavior without any collision. Purpose "collision X" requests a collision with obstacle X (X ∈ {car, pedestrian}). Purpose "lidar" requests the LiDAR to perceive an obstacle right next to the car.

Table II gives statistics about the test suite extraction for various pairs of a model and a test purpose. For each pair, columns three to six give statistics about the computation of the complete test graph: its size, generation time, and peak memory requirement. Columns seven to eleven give statistics about the generated test suite: number of test cases/scenarios, total size of the test suite (sum of states and transitions), generation time, and peak memory requirement. Column twelve gives the standard deviation (σ) for the number of test case states, showing that the test cases are of comparable size. Column 13 gives the total number of lines of the generated behavior trees.

In our experiment, we generated 9700 scenarios in total for testing AV's behavior in ten different configurations involving four TPs. Importantly, this did not require any prior knowledge about the configurations that lead to realistic scenarios. Instead, such configurations were found iteratively. Indeed, after running the behavior trees in CARLA, we observed unrealistic behaviors (e.g., actors performing discrete moves rather than continuous ones), which allowed us to improve the model by introducing a scheduler.

Another improvement concerns the translation of the discrete (grid-based) moves to the continuous environment of CARLA. When improving our model (for example by adding lane changing and turns at intersection) we had to smooth the trajectories (the motion control of the actors) for more realism. To do so, we represented in the scenario runner the sequence of grid-based moves as a sequence of waypoints the actors must pass through, following a realistic trajectory obtained by applying PID (Proportional Integral Derivative) controllers at each simulation step.

V. CONCLUSION

We proposed an approach that exploits behavioral conformance testing to automatically generate, from a formal model and test purposes, scenarios in the form of behavior trees for the Our approach makes possible the formalization (as test purposes) of safety requirements defined in standards, such as ISO 26262, and the generation of test cases, similarly to [START_REF] Makartetskiy | User-friendly) formal requirements verification in the context of ISO26262[END_REF]. The test cases can be used both for producing simulation scenarios (as we illustrated here) and for assessing the correct functioning of AVs in their environment.

Concerning future work, we plan to extend our approach to take advantage of new AV simulation standards gaining popularity, such as OpenDRIVE 3 for map representation and OpenSCENARIO 4 for scenario description (both already used by CARLA). We also plan to enhance the formal model by including near miss thresholds [START_REF] Pierson | Learning risk level set parameters from data sets for safer driving[END_REF] and the traffic-sign rules [START_REF] Bagschik | Ontology based scene creation for the development of automated vehicles[END_REF], which would make possible a refined scenario selection according to different levels of criticality.

Fig. 1 :

 1 Fig.1: Overview of the proposed approach to generate behavior trees from a configuration and a test purpose

Fig. 3 :

 3 Fig. 3: Example test purpose "collision pedestrian"

Fig. 5 :

 5 Fig.5: Excerpt of the generated behavior tree for Figure4

TABLE I :

 I Formal models of the considered configurationsTable

TABLE II :

 II Test suite extraction statistics CARLA simulator. Besides being to a large extent automatic, the generated scenarios are guaranteed to focus on specified test purposes, which provide a versatile manner of targeting corner-case situations. The approach has been illustrated on several examples of configurations (scenes and moving actors).

See also https://www.tesladeaths.com/ for a list of fatal accidents.

https://github.com/carla-simulator/scenario runner