Mickaël Alexandra Würth
email: alexandra.wuerth@inria.fr

Mickaël Binois
email: mickael.binois@inria.fr

Paola Goatin
email: paola.goatin@inria.fr

Traffic prediction by combining macroscopic models and Gaussian processes

Keywords: Macroscopic traffic flow models, Godunov scheme, parameter calibration, Gaussian process modeling, loop detector and trajectory data, travel time prediction

We propose a physics informed statistical framework for traffic travel time prediction. On one side, the discrepancy of the considered mathematical model is represented by a Gaussian process. On the other side, the traffic simulator is fed with boundary data predicted by a Gaussian process, forced to satisfy the mathematical equations at virtual points, resulting in a multi-objective optimization problem. This combined approach has the merit to address the shortcomings of the purely model-driven or data-driven approaches, while leveraging their respective advantages. Indeed, models are based on physical laws, but cannot capture all the complexity of real phenomena. On the other hand, pure statistical outputs can violate basic characteristic dynamics. We validate our approach on both synthetic and real world data, showing that it delivers more reliable results compared to other methods.

I. INTRODUCTION

Macroscopic traffic flow models are employed since several decades for state reconstruction and prediction. They consist in partial differential equations (PDEs), whose solutions can be computed by numerical schemes. Their implementation requires information provided by real data, which are typically measured by magnetic loop detectors at fixed locations. Real data are also necessary to identify model parameters. However, mathematical models may fail in capturing involved traffic situations such as congestion, lane closures or accidents, even if the true values of the calibration parameters are known [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF]. Aiming to address these shortcomings, in [START_REF] Würth | Data-driven uncertainty quantification in macroscopic traffic flow models[END_REF] we adopt the statistical framework proposed in [START_REF] Higdon | Combining field data and computer simulations for calibration and prediction[END_REF], [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] by introducing a bias term to better account for possible discrepancies between the mathematical model and reality. Following [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], we model the bias by a Gaussian process (GP), which is a classical choice when dealing with computer simulations, since it provides a flexible non-parametric framework. Once calibrated, the model can be used for traffic state estimation and prediction. The former consists in reconstructing traffic states for already realized traffic scenarios, whereas the latter deals with the prediction of the unknown future to be used for real time traffic management. In general, we can distinguish between model-driven and data-driven approaches. In the first case, physical knowledge and therefore also the calibrated parameters are used to estimate and predict the traffic state. This purely model-driven approach is often criticized to be an over-simplification of the reality, subject to model choice and calibration limitations [START_REF] Briani | Inverting the fundamental diagram and forecasting boundary conditions: How machine learning can improve macroscopic models for traffic flow[END_REF], [START_REF] Shi | A physics-informed deep learning paradigm for traffic state and fundamental diagram estimation[END_REF]. In contrast, data-driven approaches estimate the traffic states from real (historical) data, using statistical or Machine Learning methods. This alternative can deal with irregularities such as noisy data or individual driving behaviors [START_REF] Briani | Inverting the fundamental diagram and forecasting boundary conditions: How machine learning can improve macroscopic models for traffic flow[END_REF]. However, it requires a large amount of data and fails to predict nonrecurrent traffic situations, such as accidents [START_REF] Allström | Hybrid approach for short-term traffic state and travel time prediction on highways[END_REF]. Thus, a natural idea is to consider hybrid approaches for vehicular traffic determination. Currently, the so called physics-informed neural networks (PINNs) are gaining more and more attention in the literature (see e.g. [START_REF] Raissi | Physics informed deep learning (part I): Data-driven solutions of nonlinear partial differential equations[END_REF]): in PINNs, neural networks are trained to solve PDEs, whose residual is integrated in the training loss function. Analogously to our work, [START_REF] Shi | A physics-informed deep learning paradigm for traffic state and fundamental diagram estimation[END_REF] focuses on the reconstruction of vehicular traffic dynamics. However, the authors deal only with traffic estimation based on density measurements and do not consider the prediction part. An extension to flow data is considered in [START_REF] Goatin | A PINN approach for traffic state estimation and model calibration based on loop detector flow data[END_REF].

A. Contribution

In this work, we propose a hybrid method for vehicular traffic prediction, which combines the physics with GPs, as in [START_REF] Chen | APIK: Active physicsinformed kriging model with partial differential equations[END_REF]- [START_REF] Long | AutoIP: A united framework to integrate physics into Gaussian processes[END_REF]. More precisely, the observed data are modeled by a GP, which is computationally very efficient, especially when dealing with large amount of data. Moreover, we force the GP to satisfy the model PDE at virtual, i.e. unobserved points. The GP thus calibrated provides predictions for (boundary) loop detector data, which are then given to the simulator to reconstruct the space-time evolution of the traffic speed in order to compute travel times. The reconstructed traffic speeds and travel times are then compared to the ground truth. Due to limited access to both trajectory and average loop detector data, the analysis is not only performed on real world traffic scenarios but also on synthetic data generated by a microscopic simulator.

To the best of our knowledge, this hybrid approach represents an original contribution, since previous works either deal with average loop detector data prediction (see e.g. [START_REF] Mohammed | A machine learning approach to shortterm traffic flow prediction: A case study of interstate 64 in Missouri[END_REF], [START_REF] Polson | Bayesian analysis of traffic flow on interstate I-55: The LWR model[END_REF]) or travel time prediction based on purely data-driven methods (see e.g. [START_REF] Kwon | Day-to-day travel-time trends and travel-time prediction from loop-detector data[END_REF], [START_REF] Rice | A simple and effective method for predicting travel times on freeways[END_REF]).

B. Outline

The paper is organized as follows. Section II details the mathematical model and its discretization. The statistical framework for bias correction and traffic data prediction is detailed in Section III. Section IV describes the data sets considered for the validation tests, which are presented in Section V. Conclusions and perspectives are commented in Section VI.

II. DESCRIPTION OF THE MATHEMATICAL MODEL

We model the traffic dynamics by a macroscopic model consisting in a PDE that guarantees the conservation of the number of cars on the road and describes the spatiotemporal evolution of vehicle density. We refer to the Lighthill-Whitham-Richards (LWR) [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF], [START_REF] Richards | Shock waves on the highway[END_REF] model, whose key assumption, besides mass conservation, is that the average speed v is a function of the density ρ, i.e. v = V(ρ), which is referred to as the fundamental diagram. Thus, by the so called hydrodynamic relation q = ρv, the traffic flow q is also a function of the density: q = Q(ρ) = ρV(ρ). The model then consists in the scalar conservation law

∂ t ρ + ∂ x (ρV(ρ)) = 0, (1)
which is also referred to as first order model. In general, the LWR model allows to distinguish between free flow and congested traffic regimes, but it is less suitable for describing more complex situations, such as capacity drops and stop-andgo waves. Additionally, a single fundamental curve is not able to capture complex dynamics observed in congested regimes.

Several models were proposed to address these shortcomings, among which the Aw-Rascle-Zhang (ARZ) model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF], [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF] and the Generic Second Order traffic flow Model (GSOM) [START_REF] Mammar | Generic second order traffic flow modelling[END_REF]. In these models, consisiting of a system of two PDEs, the speed function depends not only on the density but also on a Lagrangian vehicle property w, which can be interpreted as an empty road velocity, measuring the more or less aggressive behaviour of the drivers. This results in a family of fundamental curves, which can capture better the spread of the data in the congested region. However, it turned out that the increase of the model dimension by the less interpretable and non-measurable variable w does not necessarily lead to better performances when applied to real traffic scenarios (see e.g. [START_REF] Shi | A physics-informed deep learning paradigm for traffic state and fundamental diagram estimation[END_REF], [START_REF] Piccoli | Second-order models and traffic data from mobile sensors[END_REF]- [START_REF] Würth | Road traffic flow reconstruction and prediction with macroscopic models enhanced by data-based statistical approaches[END_REF]). Therefore, in this work we restrict our analysis to the classical first order LWR model. Since our aim is traffic dynamics reconstruction and prediction on a road stretch, we consider the initial boundary value problem (IBVP) for (1) on a bounded interval ]x in , x out [ ⊂ R, where the variable ρ = ρ(t, x) ∈ [0, R] is equipped with prescribed initial and boundary data at t = 0 and x = x in , x = x out .

A. Numerical solution

To compute numerically approximate solutions of (1), we use the following scheme. Given a (possibly non-uniform) spatial discretization {x 0 , . . . , x M } of the interval ]x in , x out [ with x 0 = x in and x M = x out , we set the cell sizes ∆x j := x j -x j-1 for j ∈ {1, . . . , M } and a time step ∆t satisfying the Courant-Friedrichs-Lewy (CFL) condition:

∆t • max ρ∈[0,R] Q ′ (ρ) ≤ min j∈{1,...,M } ∆x j . (2)
We construct a finite volume [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] approximate solution of (1) of the form

ρ ∆xj (t, x) = ρ n j for (t, x) ∈ C n j = [t n , t n+1 [ ×[x j-1 ,
x j [ and n ∈ N. Since we also want to integrate ramps, which are modeled as junctions (see e.g. [START_REF] Garavello | Models for vehicular traffic on networks[END_REF]) we define by r n j (resp. s n j ) the measured on-ramp (resp. offramp) fluxes at position x j and time n∆t, which leads us to the formulation of the extended discrete LWR equations: if r n j ≥ 0 and s n j = 0 (and r n j-1 = s n j-1 = 0):

ρ n+1 j = ρ n j + ∆t ∆x j F n j-1 -min D(ρ n j ), max{P j S(ρ n j+1 ), S(ρ n j+1 ) -r n j } , ρ n+1 j+1 = ρ n j+1 + ∆t ∆x j+1 min{D(ρ n j ) + r n j , S(ρ n j+1 )} -F n j+1 ;
if s n j > 0 and r n j = 0 (and r n j-1 = s n j-1 = 0):

ρ n+1 j = ρ n j + ∆t ∆x j F n j-1 -min{D(ρ n j ), s n j } -min max D(ρ n j ) -s n j , 0 , S(ρ n j+1 ) , ρ n+1 j+1 = ρ n j+1 - ∆t ∆x j+1 F n j+1 -min max{D(ρ n j ) -s n j , 0}, S(ρ n j+1 ) ,
for j ∈ {2, . . . , M -1}. Above, the priority parameter P j ∈ [0, 1] is approximated by the number of lanes of cell j divided by the number of lanes of cell j added to the number of lanes of the corresponding on-ramp. Moreover, we choose the discretization in such a way that we cannot have two ramps at subsequent cell interfaces. The above scheme is based on the widely-used Godunov scheme [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] in its Cell Transmission Model (CTM) version [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF], where the fluxes across cell interfaces are given by the minimum of the sending capacity (demand D) of the upstream cell and the receiving capacity (supply S) of the downstream one. Thus, the flux F n j is computed for j ∈ {1, . . . , M -1} as

F n j = min D(ρ n j ), S(ρ n j+1 ) ,
where D(ρ) = Q(min{ρ, ρ cr }), S(ρ) = Q(max{ρ, ρ cr }) and the critical density ρ cr is given by ρ cr = argmax ρ Q(ρ).

For j ∈ {1, M }, we consider two different implementations of the boundary conditions:

1) Flow boundary conditions:

ρ n+1 1 = ρ n 1 - ∆t ∆x 1 F n 1 -min{q n in , S(ρ n 1 )} , ρ n+1 M = ρ n M - ∆t ∆x M min {D(ρ n M ), q n out } -F n M -1 ,
where q n in (resp. q n out ) denotes the inflow (resp. outflow) measured by the left (resp. right) boundary detector.

2) Density boundary conditions:

ρ n+1 1 = ρ n in , ρ n+1 M = ρ n out ,
ρ n in , ρ n out being the measured or reconstructed densities at time t = t n .

The choice typically depends on the application: for traffic flow reconstruction, the flow data usually lead to better performances. However, for travel time predictions, the density implementation results to be more favorable (see e.g. [START_REF] Würth | Road traffic flow reconstruction and prediction with macroscopic models enhanced by data-based statistical approaches[END_REF]).

In the following, we consider the Newell-Franklin [START_REF] Franklin | The structure of a traffic shock wave[END_REF], [START_REF] Newell | A theory of traffic flow in tunnels[END_REF] speed function

V(ρ) = V 1 -exp C V 1 - R ρ , (3)
whose parameters are the maximum speed V > 0, the wave propagation speed in congestion C > 0 and the maximum density R > 0, denoted by θ = (V, C, R). Note that it holds

V(ρ) ≥ 0 for ρ ∈ [0, R], V(R) = 0, V(0) = lim ρ→0 V(ρ) = V and Q ′′ (ρ) < 0 for ρ ∈ ]0, R].
Remark 1: Since in the implementation of our algorithms the measured initial and density boundary data can exceed the parameter value R, we perform a projection algorithm before executing the numerical scheme: given a density ρ > R and a speed v ≥ 0, the projected densities are computed by the inverse of the speed function (3) at v.

III. STATISTICAL METHODS FOR TRAFFIC RECONSTRUCTION AND PREDICTION

The model simulator described above can be used for reconstructing and predicting traffic scenarios by integrating the measured data y F , also called "field" observations. It is generally assumed that the field data are noisy measurements of the real quantity y P , i.e. y F (t, x) = y P (t, x) + ε at time t and position x. In the absence of any a-priori knowledge, the observation error ε is assumed to be independent and identically normally distributed (iid) with zero mean, i.e. ϵ ∼ N (0, σ 2 ε ) and σ 2 ε > 0. Kennedy and O'Hagan (KOH) [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] propose to also take into account the discrepancy between the mathematical model with optimal parameter θ * = (V * , C * , R * ) and the reality, adding a bias term b such that

y P (t, x) = y M (t, x, θ * ) + b(t, x, θ * ),
where both the simulation output y M and b depend on the parameter θ * . We note that the variables y k , k ∈ {F, P, M }, can stand for any quantity of interest, typically the flow, speed or density in the traffic context. Finally, for all i ∈ {1, . . . , N } it holds:

y F (t i , x i ) = y M (t i , x i , θ * ) + b(t i , x i , θ * ) + ε,
where X N = ((t 1 , x 1 ), . . . , (t N , x N )) denotes the set of timespace points where observations have been recorded. To estimate the bias function, we rely on a GP regression [START_REF] Higdon | Combining field data and computer simulations for calibration and prediction[END_REF], [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], which amounts to consider the discrepancy as a realization of a (zero-mean) multivariate normal distribution, i.e.

b N ∼ N (0 N , K N ) with K N = σ 2 (C N +gI N ) and g = σ 2 ε σ 2
, where b N denotes the set of observed (noisy) biases, K N (resp. C N ) the covariance (resp. correlation) matrix between these observations and l 1 , l 2 , g, σ 2 the hyper-parameters. We refer to [START_REF] Würth | Data-driven uncertainty quantification in macroscopic traffic flow models[END_REF]- [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for more details on GP modeling. Our aim is to match the real data as well as possible and to predict the system evolution in the future. We present below possible approaches for identifying θ and for forecasting traffic data.

A. Calibration

Parameter identification from measured data is a fundamental step for model validation and real world implementation. Macroscopic traffic flow models are often calibrated by fitting the fundamental diagram to data (see e.g. [START_REF] Dervisoglu | Automatic calibration of the fundamental diagram and empirical observations on capacity[END_REF]- [START_REF] Fan | A collapsed generalized Aw-Rascle-Zhang model and its model accuracy[END_REF]). However, in congested regions, traffic data are usually widely spread. Moreover, data are not necessarily measured on the whole diagram and the number of observations are often imbalanced between the free flow and congestion regimes. Thus, we focus instead on a model-driven approach which integrates the numerical solution y M in the optimization process aiming at minimizing the least squares distance among field and simulated data. The optimal parameter is then given by

θ * = argmin θ 1 N N i=1 y F (t i , x i ) -y M (t i , x i , θ) 2 . ( 4
)
Since the pure simulation output rarely fits the reality [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF], we correct the mathematical model adding the discrepancy term modeled by GPs. In particular, the bias at N new locations X N , given the observations b N , still follows a GP [START_REF] Gramacy | Surrogates: Gaussian process modeling, design, and optimization for the applied sciences[END_REF], i.e.

b( X

N) | b N ∼ N m N ( X N ), s 2 N ( X N , X N ) with m N ( X N ) = E[b( X N)|b N ] = σ 2 c N ( X N ) ⊤ K -1 N b N , s 2 N ( X N , X N ) = Cov[b( X N), b( X N)|b N ] = σ 2 c( X N , X N)-(σ 2 ) 2 c N ( X N ) ⊤ K -1 N c N ( X N
), where we choose c(•, •) as the Gaussian kernel and c N ( Xn ) = (c( X (j) N , X

N )) 1≤j≤ N ,1≤i≤N . Finally, the corrected (simulated) data y M c at time t and position x are given by y

M c ( X N , θ * ) = y M ( X N , θ * ) + m N ( X N ). (5)
Observe that the hyper-parameters necessary to compute the predictive mean are obtained by maximizing the concentrated log-likelihood function

log L(l 1 , l 2 , g, b N ) = - N 2 log 2π - N 2 log σ2 (l 1 , l 2 , g, b N ) - 1 2 log |C N + gI N | - N 2 , (6)
where the process variance is given by

σ2 (l 1 , l 2 , g, b N ) = b ⊤ N (C N + gI N ) -1 b N n .
We emphasize that this bias correction helps also to address the shortcomings of the LWR model (see e.g. [START_REF] Würth | Road traffic flow reconstruction and prediction with macroscopic models enhanced by data-based statistical approaches[END_REF]), which were mentioned in Section II.

Remark 2: Other approaches can be used for parameter identification, such as a 2-step optimization procedure of the concentrated log-likelihood function or a Markov chain Monte Carlo (MCMC) technique, as done in [START_REF] Higdon | Combining field data and computer simulations for calibration and prediction[END_REF]. However, preliminary experimental tests showed that the computationally less expensive least squares approach provides similar or even better performances after correcting the simulation by its kriging mean, as indicated in Equation ( 5).

B. Prediction

In this section, we present several approaches to predict the boundary loop detector data at future times, denoted by ŷB , which are necessary to run the numerical simulation for traffic forecast. Note that considering the first order LWR model requires the prediction of only one traffic quantity (density or flow). In general, prediction methods can either be purely data-driven, relying on historical data, or consider only data corresponding to a (short) preceding time window and which might be additionally model-driven.

1) Prediction from historical data: In the literature, we can find several prediction methods based on historical traffic data. For example, [START_REF] Kwon | Day-to-day travel-time trends and travel-time prediction from loop-detector data[END_REF], [START_REF] Rice | A simple and effective method for predicting travel times on freeways[END_REF] rely on linear regression. However, we focus on two approaches that can be considered as the most intuitive or frequently used ones. In the following, the historical traffic data form the train data set, which is scaled by first subtracting its mean and then dividing by its standard deviation. This normalization is done for each loop detector separately.

DTW: The Dynamic Time Warping (DTW) approach selects a time series in the train data set which behaves similarly to the data which immediately precede the prediction time slot. We will call these data test data, also normalized by the mean and standard deviation of the train data. To measure the similarity between two sequences, we follow [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] and use the DTW-metric DTW from MATLAB [38], which computes a predefined distance, such as the Euclidean one, between aligned time series. This alignment allows to find similarities between time series with shifted patterns or which evolve differently in time. The reference time series in the train data set is the one with the smallest DTW-distance with respect to the test data. We then perform a linear least square regression to scale and shift the reference data to the test one. The succeeding observations of the reference time series, adapted by the regression parameter, form the predicted (boundary loop detector) data.

LSTM: The long short-term memory (LSTM) recurrent neural network (RNN) is the most frequently used neural network in the context of time series analysis [START_REF] Hochreiter | Long short-term memory[END_REF]. In particular, it is capable to detect long-term dependencies between time series. To create our LSTM regression network, we make the following specifications1 : the size of the sequence input layer coincides with the number of considered loop detectors; the number of hidden units of the LSTM layer is set to 128; for the output time series, we consider a fully connected layer with the same size as the input layer and finally we include a regression layer; in the training, we use the adam-optimizer with 400 epochs and a learning rate of 0.001; the SequencePaddingDirection (resp. Shuffle) specification is set to left (resp. every-epoch). Once the network is trained on the above defined architecture, we use it to predict future time steps iteratively, transmitting nonupdated predicted values to the predictAndUpdateState function.

2) Prediction without historical data: If historical data are not available, we consider the following options to predict the future boundary data ŷB based on a (short) preceding time window. The simplest approach is to keep the data constantly equal to the last recorded measurement or an average of the last observed data, as proposed in [START_REF] Briani | Inverting the fundamental diagram and forecasting boundary conditions: How machine learning can improve macroscopic models for traffic flow[END_REF] as a comparative method. However, experimental tests show that this simplistic approach, which disregards traffic evolution information, does not lead to convincing results. This motivates us to develop approaches which exploit better the traffic dynamics of the available past data.

Pure GP: The pure GP approach belongs to the class of data-driven approaches, where traffic data are modeld by a GP. In formulas, this reads

y F (X N ) ∼ N (ȳ N , K N ) with K N = σ 2 (C N + gI N ),
where the mean ȳN is computed by taking the average of all the observed data. The covariance hyper-parameters are obtained by maximizing the concentrated likelihood function, where we replace b N by y F (X N ) in [START_REF] Shi | A physics-informed deep learning paradigm for traffic state and fundamental diagram estimation[END_REF]. The predicted data are then given by

y F ( X N ) | y F (X N ) ∼ N m y N ( X N ), (s y N ) 2 ( X N , X N ) , with m y N ( X N ) = ȳ N + k N ( X N ) ⊤ K -1 N y F (X N ) -ȳN , (s y N ) 2 ( X N , X N ) := k( X N , X N ) -k N ( X N ) ⊤ K -1 N k N ( X N ),
where k(•, •) = σ 2 c(•, •) and the constant entries of the Ndimensional vector ȳ N coincide with the ones in ȳN . Then, denoting by X NB the set of observation points at boundary loop detector positions x in and x out in the future time slot, the boundary data ŷB are given by ŷB = m y N ( X NB ). Remark 3: The choice of the mean ȳN is not evident. Since predicted data reverts typically to its prior mean [START_REF] Lee | Single nugget kriging[END_REF], we believe that ȳN is a reasonable and especially simpleto-implement choice. A more advanced suggestion can be found in [START_REF] Lee | Single nugget kriging[END_REF], where they propose the so called Single Nugget Kriging method in order to reduce the influence of the prior mean on the predictions.

Hybrid GP: Aiming at improving the prediction results using the model information, we consider methods that integrate the PDE into the GP modeling. The method proposed in [START_REF] Chen | APIK: Active physicsinformed kriging model with partial differential equations[END_REF] can be applied only to non-linear PDEs whose nonlinear term consists in products of derivatives, which is not our case. A more general approach, without restrictions on the form of the PDE, is suggested in [START_REF] Long | AutoIP: A united framework to integrate physics into Gaussian processes[END_REF]. The idea is to construct two likelihoods, a data and a virtual one. The first likelihood serves to fit the observations and the second one to fulfil the PDE equation at so called virtual points. Since there is no closed form for the posterior distribution available, they end up with a variational posterior expression. This results in solving a high dimensional optimization problem, where the number of parameters depends on the number of observations and virtual points, making the approach unsuitable for real world scenarios. This motivates us to propose a new hybrid approach, which applies to all kind of differential equations and whose set of hyper-parameters does not increase compared to the pure GP modeling. The method is based on a multi-objective optimization (MOO) of two cost functions. On one hand, we model the data by a GP, resulting in the minimization of the negative concentrated log-likelihood function:

min l1,l2,g f obj 1 (l 1 , l 2 , g) = min l1,l2,g -log L l 1 , l 2 , g, y F (X N ) .
On the other hand, we also require [START_REF] Tuo | Efficient calibration for imperfect computer models[END_REF] to be satisfied at virtual points, denoted by X Ñ = ( t1 , x1 ), . . . , ( t Ñ , x Ñ ) . This leads to the formulation of the second objective, namely the minimization of the PDE residuals at X Ñ :

min l1,l2,g f obj 2 (l 1 , l 2 , g) = min l1,l2,g ∂ t y F ( X Ñ ) + ∂ x y F ( X Ñ )V y F ( X Ñ ) ,
where y F ( X Ñ ) = m y N ( X Ñ ). The derivative expressions are computed by deriving the kernels, thus it holds, for z ∈ {t, x},

∂ z y F ( X Ñ ) = ∂ z k N ( X Ñ ) ⊤ K -1 N y F (X N ) -ȳN , and ∂ z k N ( X Ñ ) = Cov b(t i , x i ), d dz b( tj , xj )
1≤i≤N,1≤j≤ Ñ . We observe that the second objective is expressed in terms of the density. Consequently, it is natural to implement density boundary conditions in the numerical scheme. This is why we consider y F (X N ) = ρ(t 1 , x 1 ), . . . , ρ(t N , x N ) in the MOO approach (and also in the pure GP approach).

Remark 4: The choice of virtual points X Ñ is not evident and has a strong influence on the prediction. In our application, we generate uniformly distributed random points, which seems to deliver reasonable results. However, improvements might probably be achieved by considering more involved methods, as the "active PDE-informed Kriging" (APIK) approach proposed in [START_REF] Chen | APIK: Active physicsinformed kriging model with partial differential equations[END_REF].

For the multi-objective optimization purpose, we compute 100 points on the Pareto front [START_REF] Miettinen | Nonlinear multiobjective optimization[END_REF], using the MATLAB function paretosearch. We then rely on the simple kneepoint method to select a solution on the Pareto front without any prior knowledge [START_REF] Setoguchi | A knee-based EMO algorithm with an efficient method to update mobile reference points[END_REF]. The knee-point maximizes the Euclidean distance to the segment connecting the extreme points of the Pareto front, see Figure 1. It is considered as a reasonable solution since moving along the Pareto front would lead to a larger deterioration in one of the objectives. Once the knee-point is determined, the optimal hyperparameters l 1 , l 2 and g are identified. This enables us to compute the desired boundary loop detector data by exploiting the predictive mean formula, namely ŷB = m y N ( X NB ). Remark 5: The second objective function f obj 2 reminds of the residual loss function in the PINNs approach [START_REF] Raissi | Physics informed deep learning (part I): Data-driven solutions of nonlinear partial differential equations[END_REF], where the PDE solution is approximated by a neural network trained on a large amount of so called auxiliary points. In our case, the GP predictive equations are computationally less expensive. Moreover, we keep the two cost functions separated, instead of including them in a single objective [START_REF] Long | AutoIP: A united framework to integrate physics into Gaussian processes[END_REF].

C. Determination of travel times

Aiming to predict travel times in some "future" time interval [T now , T F ], we integrate the previously presented approaches in the following procedure:

1) Calibration of θ by the least squares approach in some preceding time window [T I , T now [.

2) Estimation of predicted boundary data ŷB in

[T now , T F ].
3) Computation of of the model solution y M and eventually its bias correction b in

[T now , T F ]× ]x in , x out [. 4) Estimation of travel times in [T now , T F ].
In particular, the output of the simulation is used to calculate travel times by solving the ordinary differential equation

dx dt = v(t, x(t)), x(t 0 ) = x in , for t 0 ∈ [T now , T F ].
The travel time is then given by the first time τ = τ (t 0 ) > t 0 such that x(τ ) = x out . The whole procedure is summarized in Algorithm 1.

Algorithm 1 Travel time computation from simulation output. Require: Departure time t 0 of a vehicle starting at position

x in and going to position x out , simulated speed v(t, x) at time t and position x, time step size ∆t of numerical simulation. Initialize space position by x = x in ; Initialize travel time by τ = 0; while x < x out do Update space position by x = x + ∆t • v(t 0 + τ , x); Update travel time by τ = τ + ∆t; end while Return travel time τ .

To evaluate the performance of our travel time prediction, we need a comparison to real data, which in the best case are recovered from Global Positioning System (GPS) tacking, probe vehicles or video recordings. However, these measurements are not often available or accessible [START_REF] Kwon | Day-to-day travel-time trends and travel-time prediction from loop-detector data[END_REF]. In these cases, we will use aggregated loop detector measurements as an approximation to the ground truth. As an error metric, we consider a relative version of the root mean square error (rRMSE), comparing the true and reconstructed travel times, τ and τ , for the same road stretch at N τ different departure instants. Thus, the total travel time error E τ reads as

E τ = Nτ i=1 (τ i -τi ) 2 Nτ i=1 τ 2
i . In the following, we detail the computation of the reference travel times in the two mentioned cases: a) Trajectory data: Real travel time data may fluctuate a lot due to different traffic situations such as congestion and free flow regimes, different driving behaviors and vehicle types or due to external factors such as weather conditions. Instead, the mathematical model provides average speed values and therefore only average travel times can be obtained by Algorithm 1. For a better comparison, we also average the real recorded data: to compute the reference travel time for a vehicle starting at time t, we take the mean over the travel times of vehicles which started their trip between t -ε and t + ε, where ε > 0 has to be chosen.

b) Aggregated data: If no or only few trajectory data are available, it is possible to derive travel times approximations from aggregated data using the N-curves method [START_REF] Treiber | Traffic flow dynamics: data, models and simulation[END_REF], obtained by summing up the aggregated flow data measured by loop detectors. To estimate the travel time between two loop detectors A and B for a vehicle starting at time t A at loop A, we calculate the number of accumulated vehicle counts for loop A at time t A and we intersect this number with the N-curve of loop B, obtaining the arrival time t B . The travel time is thus given by τ = t B -t A . Due to the usage of aggregated data, this approach leads directly to average travel times. However, a major source of error for the N-curve method is the presence of on-of off-ramps, lane changes and overtaking maneuvers [START_REF] Treiber | Traffic flow dynamics: data, models and simulation[END_REF]. An alternative for travel times reconstruction can be the application of Algorithm 1, where the entries of the field v(t, x) are given by the piece-wise constant average speed measurements. We refer to this as the baseline method.

IV. DESCRIPTION OF TRAFFIC DATA SETS

To validate the proposed approaches, we refer to traffic data recorded by loop detectors, providing aggregated information over time, such as the traffic flow q and the occupancy O (the percentage of time a detector is occupied by a vehicle [START_REF] Polson | Bayesian analysis of traffic flow on interstate I-55: The LWR model[END_REF]). From the occupancy measurement, one can derive the traffic density ρ from the relation ρ = O/l, where l denotes the average vehicle length. If two detectors are installed in very close succession, the average traffic speed v can be also directly computed. Otherwise, it can be derived from the fundamental relation v = q/ρ, which leads to a spatial and not a temporal average value. We point out that contrary to the flow (resp. density) which belongs to the class of temporal (resp. spatial) traffic data, the speed can be defined as both a temporal and a spatial quantity [START_REF] Treiber | Traffic flow dynamics: data, models and simulation[END_REF]. However, these two definitions differ from each other, thus they naturally lead to different results in applications such as travel time predictions.

A. Synthetic microscopic traffic data

Since real data often suffer serious limitations, such as nonfunctioning sensors or measurement errors, we first validate the proposed approaches on synthetic data. To this aim, we rely on data generated by Simulation of Urban MObility (in short SUMO) [START_REF] Lopez | Microscopic traffic simulation using SUMO[END_REF]. Unlike data derived by simulations based on a macroscopic traffic flow model [START_REF] Würth | Validation of calibration strategies for macroscopic traffic flow models on synthetic data[END_REF], the microscopic simulator tracks each vehicle, giving direct access to trajectory data and the corresponding travel times. We consider two traffic scenarios, both simulating a highway traffic situation for a rampless 10km road stretch with three lanes and a constant speed limit of 100km/h, equipped with 3 × 10 loop detectors, one for each lane at ten different, nonequidistant locations. Moreover, the sensor data are aggregated over 6 minutes. The traffic flow consists of three different vehicle types differing in their desired maximum speed. The length of all vehicles is set to 5m. SUMO provides all measurable traffic quantities: the flow, occupancy and speed per lane. To compute the average speed over the three lanes, we use the space mean formula v = 1 ρ 3 i=1 ρ i v i , v i given by the fundamental equation. We note that other choices for v i are possible, such as the arithmetic or harmonic mean speed; however, for these alternatives, we observed worse traffic speed reconstructions.

SUMO-1: The first synthetic scenario reproduces a highly congested traffic situation, generated by gradually reducing the speed limit on the last 0.01km of the road, thus inducing a backward moving congestion wave, which dissipates after the considered time window of 3 hours. In Figure 2a, we depict the speed values measured by the loop detectors for the 6 minute aggregated time intervals. The distribution of all the data points is visualized by the fundamental diagrams in Figure 2b. The non-monotonicity of the flow data indicates the presence of different traffic regimes: the free flow and congestion. In particular, the congested phase presents more widely spread data than the free flow, which is typical of real traffic scenarios. For computing reference travel times, we refer to Section III-C a), where we choose ε = 10s . The first (resp. last) considered departure time is after 6 minutes (resp. 2 hours and 30 minutes), resulting in N τ = 865 vehicle trajectories. Considering Figure 5a, we observe a steadily increasing travel time up to τ ≃ 27min due to congestion. The results for the method of N-curves (resp. baseline method) is depicted in black (resp. blue). It appears that the N-curve method reflects better the ground truth than the baseline one, which is confirmed by the rRMSEs: E τ N-curve = 0.027 < E τ baseline = 0.052. This can be explained by the more precise information on flow variations exploited by the N-curves, compared to the piece-wise constant speed data employed by the baseline approach.

SUMO-2: The second scenario simulates an accident occurring at time t = 1h. We implement this in SUMO by closing the rightmost lane between kilometer 6 and 9 for a duration of 80 minutes. This creates less congestion than in SUMO-1, see Figure 3a and the few data points in the congested region in Figure 3b. However, we observe a sudden change in the traffic regime, especially for the loop detector located at x = 5.5km, right before the lane closure. This scenario is interesting to analyze since the sudden change in the dynamics cannot be captured by the model without ad-hoc changes. Moreover, the last boundary detector is not affected by this lane closure. In this case, the role of the bias becomes crucial.

B. Real traffic data

Aiming to test our approach in real world scenarios, we consider the RTMC data set [START_REF]Mn/Dot Traffic Data[END_REF], provided by the Minnesota Department of Transportation (MnDOT). As the synthetic SUMO data, the RTMC data are 6 minute averages obtained by single loop detectors measuring the traffic flow and the occupancy. For the tests, we consider a 4.85km long road stretch on the northbound direction of the interstate highway I-35W, which is equipped with 8 sensors on the mainlane, 2 at on-ramps and 3 at off-ramps. The road stretch counts five lanes and the speed limit is 55miles/hour (≈ 90km/h). We refer to the year 2013 measures, serving as the historical train data, and some selected days in 2014 as test data. The data were pre-processed to remove abnormalities: we discarded eight full days (due to the presence of negative traffic quantities) and we replaced all measured flow (resp. density) values by zero if the corresponding density (resp. flow) was zero. In these cases, we set the average speed to 200km/h. We remark that a zero flow could also correspond to a fully congested road, however this is never the case for our extracted RTMC data, which can be easily verified by comparing the data measured by nearby loop detectors.

RTMC: The selected real data scenario covers the morning time slot from 6am to 9am of Wednesday, November 5th, 2014. Figure 4a show free flow conditions in the first 30 minutes. Then, the traffic gets denser until almost the end of the considered time period, see also the flow fundamental diagram in Figure 4b. Again, the congested part show a larger spread of the data for densities higher than 110veh/km. Since in this case we do not have access to trajectory data, the reference travel times can be only estimated by using the aggregated loop detector data. We remark that the initialization of the N-curves belonging to the first and last loop detector are based on a free flow assumption at 6am. In Figure 5b, we illustrate the results of the two approaches for vehicles departing during the first 2 hours. Unlike the SUMO scenario, here the N-curve method fails since the travel time decreases constantly, resulting finally in negative values. The baseline approach however captures the increasing of travel times, due to lower measured speeds in the second part of the time slot.

V. MODEL AND STATISTICAL METHODS VALIDATION

Since we are interested int travel time prediction, we consider the speed as our quantity of interest. In the numerical scheme, we focus on a density boundary implementation, which turned out to perform better in our context. Moreover, the initial datum is approximated by piece-wise constant density values, which are given by the first recorded aggregated measurements.

A. Traffic calibration and reconstruction

The calibration is run on the first 2 out of 3 hours for each scenario, using the MATLAB nonlinear optimization function fmincon. The last hour is left for the prediction tests in Section V-B. In Table I, we list the parameters and the rRMSE E (resp. E c ) at observation points between the field and simulated (resp. corrected) speed data. For comparison, and to underline the benefit of the physical model, we consider the error E GP between the kriging mean of y F (X N ) ∼ N (ȳ N , K N ) and the field data. We observe that the bias corrected version E c always outperforms E and E GP . Moreover, the errors for the We recall that the RTMC scenario includes ramps; however, we set r n j = 0 = s n j in the simulation code, since this choice leads to slightly lower rRMSEs, compared to E = 0.228, E c = 0.071 and θ * = (107, 86, 290) including ramp information. We will proceed analogously in the prediction section, thus avoiding to forecast ramp data, whose impact is questionable.

We also point out the strong improvement of the error metric after the bias correction in SUMO-2. As previously discussed, the pure simulation cannot reflect the traffic dynamics properly, because the model is not designed to capture the lane closure. This scenario emphasizes the benefit of the bias modeling, which is also underlined by the space-time speed plots in Figure 6 (a) and (b). Indeed, the corrected version is able to detect the traffic jam. This can also be observed in the speed profile for the sixth sensor, see Figure 6 (c).

B. Traffic prediction

For the prediction of the density boundary data needed to run the numerical scheme, we introduce the density boundary prediction rRMSE Êρ B between the measured aggregated y F ρ ( X NB ) and predicted ŷB density data at 2n t future boundary observation points. We also compare the simulated (resp. corrected) speed with the available coarse field data at loop detector positions and aggregated future times by the rRMSE, denoted by Ê (resp. Êc ). These metrics are used as indicators for travel time prediction quality, lower values corresponding to better forecasts. Additionally, we compute the two speed metrics using the real boundary data y F ρ ( X NB ) in the simulation. This will be referred to as the oracle case. Finally, in the SUMO data case, we also compare the predicted relative travel time error metrics Êτ and Êτ c with the reference data, where we consider the vehicles departing between t = 90min and t = 150min, corresponding to 361 trajectories.

Remark 7: For the determination of the number of virtual time points which are used in the second objective function, we differentiate between the time and space dimension: for the first one, we double the size of the 6 minute aggregated data in the 3 hours, thus sampling 180 6

• 2 = 60 points from the uniform distribution. The virtual space points however are not sampled, instead they match exactly the position of the main loop detectors. We observe that sampling the space points lead to worse results in our case.

In Table II, we highlight in bold the lowest value for each considered error metric. First of all, we observe that, compared to the previous estimation section, the difference between the predicted speed errors Ê and Êc are negligible, meaning that the bias correction in the future time slot has almost no or even a negative impact. Second, a good density boundary prediction does not necessarily lead to the lowest speed error Ê, as in the oracle boundary case where it holds Êρ B = 0 but Ê can be bigger than in other approaches. This underlines the difficulty when dealing with different traffic quantities: since the model describes density dynamics, a good reconstruction of the density does not necessarily result in a good speed estimation. In this regard, the independent speed reconstruction provided by second order models do not lead generally to improvements, as shown in [START_REF] Würth | Road traffic flow reconstruction and prediction with macroscopic models enhanced by data-based statistical approaches[END_REF]. Third, against expectations, the lowest Ê value does not correspond to the best travel time prediction result: in SUMO-1 the MOO approach outperforms the other cases although its coarse speed reconstruction is worse. In SUMO-2 all proposed methods perform very in terms of travel time prediction: the oracle boundary case shows a slight outperformance, although its speed rRMSEs are worse than the others. Looking at the traffic volume in the prediction hour in Figure 3a, the good results are easily explainable: the congestion induced by the lane closure has almost no impact on the last hour, consequently it is easier to predict the boundary data. For RTMC, we observe that the DTW approach, which provides clearly the lowest Êρ B , performs the worst in terms of speed rRMSE. In contrast, although the boundary predictions obtained by the LSTM network are disappointing, its speed prediction power is convincing.

In Figure 7, we illustrate the results for SUMO-1 and the MOO approach. Although the knee-point heuristic does not lead to the lowest possible density boundary prediction error Êρ B (see Figure 7a), the performance is acceptable since only a few points on the Pareto front undercut this point. Moreover, in order to understand the rather disappointing performance of the oracle case, we illustrate in Figure 7b the travel time profiles for the simulations and oracle cases. By the black vertical line, we highlight the starting point of the 361 considered trajectories, used in the travel time error computation. Naturally, the oracle boundary curve and corrected simulated curve coincide in past because the is executed with the same the last hour, the oracle boundary case starts to overestimate the travel times. In contrast, the (corrected) simulation both underand overestimates the speed which leads therefore to a lower rRMSE in average. Due to the increasing uncertainty in the future, it is clear that the correction tends to go back to the pure simulation, which underestimates the speeds when using real boundary data. This in turn explains the behavior of the oracle curve and the similar performances of Ê and Êc .

As a final analysis, we compare in Figure 8 the uncertainty of the pure GP and our proposed MOO approach. Together with the real future boundary observations (red stars) and its predictions ŷB = m y N ( X NB ) (blue line), we depict 100 realizations of N m y N ( X NB ), (s y N ) 2 ( X NB , X NB ) and the corresponding 90% predictive intervals. We observe that in the hybrid MOO approach, all observations are lying inside the predictive intervals, whereas in the pure GP approach the uncertainty is higher for the right boundary detector. This last analysis reinforces the evidence of the advantage of integrating the physics in the prediction approach.

VI. CONCLUSION

In this paper, we focused on traffic flow reconstruction and prediction by using a first order macroscopic traffic flow model and statistical approaches. We have shown the benefit of correcting the simulation output by a bias in order to compensate model limitations in reconstructing real data. For prediction, bias correction can be neglected due to its (a) Pure GP approach. convergence towards the pure simulation. Moreover, we proposed a hybrid MOO approach for boundary data forecast, combining the physical knowledge, given by the mathematical model, and GPs. This method is based on multiobjective optimization and it addresses the shortcomings of the existing ones: it applies to all kind of differential equations and the set of hyperparameters does not increase compared to the pure GP modeling. The results are robust, providing competitive error metric values in all tested scenarios. Additionally, we observed a non-logical anti-correlated behavior of the metrics Ê and Êτ , which might be partly due to the choice between time or space mean speeds. Indeed, the difficulty of predicting the speed is a frequently mentioned problem (see e.g. [START_REF] Liu | A Gaussianprocess-based data-driven traffic flow model and its application in road capacity analysis[END_REF], [START_REF] Wang | Macroscopic traffic flow modelling of large-scale freeway networks with field data verification: State-of-theart review, benchmarking framework, and case studies using metanet[END_REF]). This work opens several perspectives for future research, as the choice of the virtual points in the MOO approach, for example considering the so called APIK approach [START_REF] Chen | APIK: Active physicsinformed kriging model with partial differential equations[END_REF]. Moreover, it would be interesting to compare the performance of our proposed MOO approach with the PINNs method.

Fig. 1 .

 1 Fig.1. Illustration of the knee-point method for the MOO approach.

Fig. 2 .

 2 Fig. 2. Traffic scenario SUMO-1.

Fig. 3 .

 3 Fig. 3. Traffic scenario SUMO-2.

Fig. 4 .

 4 Fig. 4. Traffic scenario RTMC.

Fig. 5 .

 5 Fig. 5. Illustration of average travel times.

Fig. 6 .

 6 Fig. 6. SUMO-2. Illustration of space-time speed and speed profile.

Fig. 7 .

 7 Fig. 7. SUMO-1. Travel time prediction results for MOO approach.

Fig. 8 .

 8 Fig. 8. SUMO-1. Uncertainty quantification for pure GP and MOO approach.

https://fr.mathworks.com/help/deeplearning/ug/time-series-forecasting-usi ng-deep-learning.html. Accessed on 06/13/2023.

Manuscript received , 2023; revised . This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.


Thomas Schatzmayr

Welp Sa

Sandrine Marceau

Clément Piegay

Fouad Laoutid

Philippe Gle

Cesar Segovia

Etienne Gourlay

Emmanuel Gourdon

Optimized vegetal wools for indoor comfort: coupling fire treatment with acoustic and hygrothermal performances

Keywords: sustainable building materials, vegetal wools, fireproof treatment, acoustics, hygrothermics 1

Vegetal wools are used as an alternative to traditional wall insulator materials, as they present great hygrothermal and acoustic properties, therefore providing interior building comfort. However, one major concern regarding the material's flammability arises, as generally, bio-based materials are poorly classified in terms of flame resistance. Therefore, fireproof treatments are needed and the knowledge regarding their influence on the wool's microstructure, and consequently, on its acoustic and hygrothermal performance, is lacking. This article presents the state of the art and the methodology that will be used to choose and optimize the fireproof treatment on the vegetal wools, and comprehend how it affects the material's microstructure and its acoustic and hygrothermal performances. In order to optimize the performance of the vegetal wools, the final objective of this project will be to develop models to describe the correlation between their micro and macro properties.

INTRODUCTION

The Intergovernmental Panel on Climate Change (IPCC) reported that the global warming effects are happening faster each year, and in order to change this scenario, urgent actions must be taken [START_REF] Gourdon | Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. For the building sector, the development of efficient and resilient buildings aiming to decarbonize the industry is a key solution as this sector is responsible for more than one third of the total carbon dioxide emitted and energy consumed in the world, per year, considering pre-pandemic data from 2019 [START_REF]Global Status Report for Buildings and Construction: Towards a Zero-emission[END_REF]. In addition to the issues of climate change and energy consumption, the noise pollution should also be addressed, as it is the second highest cause of health problems for citizens, after air pollution [START_REF] Who | Environmental noise guidelines for the European Region[END_REF]. Therefore, efficient wall insulation can contribute to the peoples' health and comfort.

In order to tackle both of these problems, bio-based insulators emerges as a candidate as they temporarily store carbon and are great acoustic and thermal insulators [START_REF] Levasseur | Assessing temporary carbon sequestration and storage projects through land use, land-use change and forestry: comparison of dynamic life cycle assessment with ton-year approaches[END_REF]- [START_REF] Arenas | Eco-Materials with Noise Reduction Properties[END_REF]. Among several possibilities of bio-based insulation materials, this work focuses on vegetal wool. This material presents interesting characteristics for a wall insulation application, such as its physical, hygrothermal and acoustic properties. However, one major concern arises regarding the material's fire resistance, especially for this application, as bio-based materials usually present low resistance against fire [START_REF] Dorez | Effect of cellulose, hemicellulose and lignin contents on pyrolysis and combustion of natural fibers[END_REF]. Therefore, to expand the applications of vegetal wools for green building envelope a fireproof treatment for this material should be considered. On the other hand, the application of a fireproof treatment on the vegetal wools could be modifying its microstructure [START_REF] Antoun | Renewable phosphorous-based flame retardant for lignocellulosic fibers[END_REF]- [START_REF] Taib | Current progress of biopolymer-based flame retardant[END_REF]. As their hygrothermal and acoustic properties are directly related to their microstructure the fireproof treatment will, consequently, affect the material's hygrothermal and acoustic performances. Thus, this should be evaluated, in order to ensure that the treated vegetal wools would still be considered a good alternative for the development of green buildings.

This paper presents the methodology that will be used to identify eco-friendly fireproof treatments and study their impact on the thermal and acoustic performances of vegetal wools. The three main objectives are:

1. Evaluate the alternatives and the application methods of the fireproof treatment in the vegetal wools; 2. Comprehend the influence and optimize the fireproof treatment on the fibers and microstructure of the vegetal wools; 3. Adapt micro-macro modeling methods to simulate the impact of the fireproof treatment on acoustic and hygrothermal properties of the vegetal wools.

In this article, the current state of the art on vegetal wools and their properties is presented, as well as the main mechanisms of flame retardancy of fireproof treatments and the main methods used to increase flame retardancy of vegetal fibers. The proposed methodology for this study is then detailed.

STATE OF THE ART

Vegetal wools

Vegetal wools are manufactured products made of lignocellulosic and synthetic fibers.

The manufacturing process can be done by three different methods: mechanical bonding (needle-punching), chemical bonding and thermal bonding [START_REF] Roy | Thermal Bonded Nonwoven -An Overview[END_REF]. Among the options, thermal bonding is the most popular for nonwovens manufacture, at it provides high production rates, while demanding less energy than chemical bonding, for example [START_REF] Roy | Thermal Bonded Nonwoven -An Overview[END_REF]. Hence, the focus of this paper will be specifically on thermal bonding process.

In this process, the vegetal fibers are mixed with thermoplastic polymeric ones, and are then heated until the melting or softening temperature of the synthetic fibers, followed by subsequent cooling to solidify and form a bond between the vegetal fibers [START_REF] Roy | Thermal Bonded Nonwoven -An Overview[END_REF]. Then, the vegetal wools obtained are bi-component materials where the thermoplastic fiber acts as a binder and the lignocellulosic fiber is the base fiber. For this study, hemp and flax will be used as the base fiber, and polylactic acid (PLA) as the binder.

The vegetal wools are generally composed by 80 to 90 wt.% of vegetal fiber, and 10 to 20 wt.% of polymeric fiber [START_REF] Piegay | Approche conjointe acoustique et thermique pour l'optimisation des laines végétales du bâtiment[END_REF]. For example, in Freivalde et al. [START_REF] Freivalde | Flammability of raw insulation materials made of hemp[END_REF] hemp wools were made of 92% of hemp fibers and 8% of polyethylene-polypropylene fibers. Fig. 1a and Fig. 1b shows the visual aspects of a vegetal wool and a SEM observation of polymer and vegetal fibers, respectively. For the manufacturing process of flexible insulation panels traditionally used in the building field, the material generally presents a low density of less than 50 kg.m -3 and a high porosity greater than 90% [START_REF] Freivalde | Flammability of raw insulation materials made of hemp[END_REF]- [START_REF] Piégay | A self-consistent approach for the acoustical modeling of vegetal wools[END_REF].

Acoustic and hygrothermal properties of vegetal wools have been measured and evaluated by a few studies, showing promising results, but further studies are needed to consolidate the knowledge around it.

For the acoustic properties research works have been done not only on vegetal fibers and shives [START_REF] Berardi | Acoustic characterization of natural fibers for sound absorption applications[END_REF]- [START_REF] Bousshine | Acoustical and thermal characterization of sustainable materials derived from vegetable, agricultural, and animal fibers[END_REF], but also on vegetal wools [START_REF] Piégay | A self-consistent approach for the acoustical modeling of vegetal wools[END_REF], [START_REF] Piégay | Acoustical model of vegetal wools including two types of fibers[END_REF]. Fig. 2a shows some results of sound absorption from hemp wool with different thicknesses. It underlines the high sound absorption potential of these materials in a wide range of frequencies, sustaining the argument that it is a good acoustic insulator.

Regarding the hygrothermal properties, some studies explored these properties in vegetal wools [START_REF] Collet | Water vapor properties of two hemp wools manufactured with different treatments[END_REF], [START_REF] Piégay | A cylindrical self-consistent modelling of vegetal wools thermal conductivity[END_REF], [START_REF] Korjenic | Development and performance evaluation of natural thermal-insulation materials composed of renewable resources[END_REF]- [START_REF] El-Sawalhi | Estimation of the thermal conductivity of hemp based insulation material from 3D tomographic images[END_REF]. On the hygric point of view, Collet et al. [START_REF] Collet | Water vapor properties of two hemp wools manufactured with different treatments[END_REF] showed that hemp wools can absorb about ten times more water than mineral wools (comparing to the literature), and that both materials show similar water vapor permeability. Fig. 2b shows the sorption curve of a hemp wool sample. For their thermal properties, hemp and flax wools shows values of thermal conductivity ranging usually between 40 and 55 mW.m -1 .K -1 [START_REF] Piégay | A cylindrical self-consistent modelling of vegetal wools thermal conductivity[END_REF], [START_REF] Korjenic | Development and performance evaluation of natural thermal-insulation materials composed of renewable resources[END_REF]- [START_REF] El-Sawalhi | Estimation of the thermal conductivity of hemp based insulation material from 3D tomographic images[END_REF], depending on the wools' density and porosity. To sum up, vegetal wools have a great potential to be used as a sustainable insulator material for green building applications, being, therefore, an alternative to traditional building materials [START_REF] Gumanová | Experimental Survey of the Sound Absorption Performance of Natural Fibres in Comparison with Conventional Insulating Materials[END_REF].

Fireproof treatments

Despite their performances adapted to the building sector, there is still a concern regarding the flammability of vegetal wools, as, indeed, bio-based materials are usually poorly classified regarding their fire resistance [START_REF] Amziane | Bio-aggregates Based Building Materials: State-of-the-Art Report of the RILEM Technical Committee 236-BBM[END_REF], [START_REF] Kumar | Comparative analysis of building insulation material properties and performance[END_REF]. Therefore, the application of a fireproof treatment could improve the vegetal wool's resistance to fire, and allow it to be applied in buildings without introducing any hazard.

Commonly used flame retardants (FR) are mineral, halogenated, phosphorousbased, silicon-based and nitrogen-based products [START_REF] Laoutid | New prospects in flame retardant polymer materials: From fundamentals to nanocomposites[END_REF]- [START_REF] Costes | Bio-based flame retardants: When nature meets fire protection[END_REF]. For vegetal fibers, great results have been achieved through grafting of phosphorous into the fibers composition [START_REF] Hajj | Grafting of phosphorus flame retardants on flax fabrics: Comparison between two routes[END_REF], [START_REF] Hajj | Influence of lignocellulosic substrate and phosphorus flame retardant type on grafting yield and flame retardancy[END_REF], [START_REF] Taibi | Flame retardancy of flax fibers by pre-irradiation grafting of a phosphonate monomer[END_REF], [START_REF] Antoun | Renewable phosphorous-based flame retardant for lignocellulosic fibers[END_REF]. Phosphorous-based FR promotes the formation of a carbonized layer, called char, which acts as a physical barrier, reducing both the thermal and oxygen diffusion, as well as the volatilization of combustible products [START_REF] Taib | Current progress of biopolymer-based flame retardant[END_REF], [START_REF] Laoutid | New prospects in flame retardant polymer materials: From fundamentals to nanocomposites[END_REF], [START_REF] Costes | Bio-based flame retardants: When nature meets fire protection[END_REF]. Furthermore, they are not toxic [START_REF] Taibi | Flame retardancy of flax fibers by pre-irradiation grafting of a phosphonate monomer[END_REF]. With bio-based materials, the use of phosphorousbased products, associated with nitrogen-based products forms intumescent flame retardant systems (IFRS). The IFRS forms an expanded char during thermal degradation, which acts as an insulating barrier, reducing heat transfer between the heat source and the material's surface, and limiting the fuel transfer and the oxygen diffusion [START_REF] Laoutid | New prospects in flame retardant polymer materials: From fundamentals to nanocomposites[END_REF]. For the formation of a IFRS, three components are required: an acid source, provided by the phosphorous thermal decomposition; a carbonizing agent, which is the bio-based material; and a blowing agent, which can be obtained from nitrogen products [START_REF] Laoutid | New prospects in flame retardant polymer materials: From fundamentals to nanocomposites[END_REF], [START_REF] Costes | Bio-based flame retardants: When nature meets fire protection[END_REF], [START_REF] Gaan | Effect of phosphorus and nitrogen on flame retardant cellulose: A study of phosphorus compounds[END_REF].

To graft the phosphorous on the fiber's molecules, different techniques might be used, in particular, chemical grafting and silane coating. The chemical grafting aims to form covalent bonds between the phosphorous and the chemical components of the fiber. If the treatment is performed in mild conditions (low temperature, short duration), a condensation reaction occurs between the phosphorous and the aromatic groups present in lignin. On the other hand, in severe conditions (high temperature, long duration), a phosphorylation reaction occurs, grafting the phosphorous on cellulose [START_REF] Antoun | Renewable phosphorous-based flame retardant for lignocellulosic fibers[END_REF], [START_REF] Hajj | Grafting of phosphorus flame retardants on flax fabrics: Comparison between two routes[END_REF], [START_REF] Hajj | Influence of lignocellulosic substrate and phosphorus flame retardant type on grafting yield and flame retardancy[END_REF], [START_REF] Moussa | Toward the cottonization of hemp fibers by steam explosion. Flame-retardant fibers[END_REF].

Alternatively, the silane treatment aims to form a coating layer on the fiber's surface, creating a thermal shield and reducing oxygen diffusion [START_REF] Castellano | Synthesis and characterization of a phosphorous/nitrogen based sol-gel coating as a novel halogen-and formaldehyde-free flame retardant finishing for cotton fabric[END_REF]. When associated with phosphorous and nitrogen products, it can form a IFRS with enhanced performance [START_REF] Decsov | Microfibrous cyclodextrin boosts flame retardancy of poly(lactic acid) II -phosphorous silane treatment further enhances the effectivity[END_REF]. Nevertheless, the work carried out with these treatments concerns applications on small quantities of fibers. It is therefore essential to evaluate whether these methods are viable for the treatment of lignocellulosic fibers on a larger scale, particularly in the context of industrial manufacturing.

Moreover, as a side-effect, the incorporation of phosphorous-based FRs can, potentially, increase the fungus development in the vegetal wools, because phosphorous is one of the possible elements that stimulates the microbial growth [START_REF] Vasiliauskienė | Impact of Cellulolytic Fungi on Biodegradation of Hemp Shives and Corn Starch-Based Composites with Different Flame-Retardants[END_REF], [START_REF] Alexopoulos | Fungus." Encyclopedia Britannica[END_REF]. This will therefore be another important aspect to consider for the flame retardant overall effectiveness.

METHODOLOGY

Aiming to develop optimized fireproof treated vegetal wools for insulation, the proposed methodology scheme for this work is presented on Fig. 3. For this study, two types of vegetal fibers will be considered: flax and hemp fiber. Both fibers are largely available in France, and have been used for flax and hemp wools production. The first step of this study will consist on understanding the manufacturing process of turning the vegetal fibers into vegetal wools as described in section II. Then, a study about fireproof treatment of the vegetal wools will be carried out, with the objective to define which flame retardant and application method that provide the desired performance. Three approaches will be evaluated: grafting with mild conditions; grafting in severe conditions; and silane treatment. The FRs considered are (3trimethoxysilylpropyl)diethylenetriamine (TRIAMO), dihydrogen ammonium phosphate (DAP), phytic acid (PA) and urea. These compounds were chosen aiming to form IFRS and to be viable for a usage in large scale. The fire properties will be evaluated through cone calorimeter test and thermogravimetric analysis. Moreover, as the main substances considered for the fireproof treatment are phosphorous-based, components that may be nutrients for fungus development, a specific study will be conducted to assess their impact on fungal growth.

In the following, it will be necessary to comprehend the physicochemical and the microstructural properties of the vegetal wools, due to the fact that the application of a fireproof treatment on the fibers could change some of its properties on both micro and meso scales. For example, fireproof treatment could impact important parameters such as porosity, airflow resistivity, density and fiber radii. Consequently, macro-scale properties such as acoustic and hygrothermal, could also be affected. This step will be fundamental to the fully comprehension and control of the vegetal wools' acoustic and hygrothermal performances. It will be done by image analysis on the optical and electronic microscopy, pycnometry and airflow resistivity meter.

Next, the acoustic and hygrothermal parameters of treated and non-treated vegetal wools will be evaluated and compared. This comparison will allow the enrichment of the existing database for untreated vegetal wools, and to provide new knowledge on the treated material. The acoustic properties can be obtained with Kundt tube characterizations and the hygrothermal properties can be obtained through the dynamic vapor sorption (DVS) test, the moisture buffer value (MBV), water vapor permeability (WVP) measurement and hot disk test [START_REF] Collet | Water vapor properties of two hemp wools manufactured with different treatments[END_REF], [START_REF] Amziane | Bio-aggregates Based Building Materials: State-of-the-Art Report of the RILEM Technical Committee 236-BBM[END_REF].

This work also includes the adaptation of modelling methods to describe the hygrothermal and acoustic behavior of fireproof treated vegetal wools. The use of these models will permit numerous simulations to be carried out by varying the input parameters and will also be useful for identifying possible improvements that could be done to optimize the treated vegetal wool performance. To do so, there are many different possible approaches, and for this study, the micro-macro approach was chosen, using the self-consistent method (SCM). This method consists on an homogenization between microscopic scale parameters and macroscopic scale properties [START_REF] Piégay | A cylindrical self-consistent modelling of vegetal wools thermal conductivity[END_REF]. Through this method, it is possible to create direct analytical relationships between micro and macro scales and to carry out rapid inverse analysis [START_REF] Piégay | A self-consistent approach for the acoustical modeling of vegetal wools[END_REF].

Piégay et al. [START_REF] Piégay | A cylindrical self-consistent modelling of vegetal wools thermal conductivity[END_REF] applied this method to develop a model for the thermal conductivity of untreated vegetal wools and achieve a mean relative deviation from experimental results of less than 2%. For the acoustics, Piégay et al. [START_REF] Piégay | A self-consistent approach for the acoustical modeling of vegetal wools[END_REF] used the SCM method to obtain the sound absorption of untreated vegetal wools. This study would incorporate the changes in the microstructure with the use of flame retardants into the modelling approach, and validate it by comparing with the previous experimental tests. Based on these modelling methods, it will be possible to manufacture optimized fireproof treated vegetal wools panels prototypes.

CONCLUSION

In this study, a novel fireproof treated vegetal wool will be developed. Studies about the physicochemical and microstructure properties, as well as the modelling methods, will be used to understand and optimize the eco-friendly fireproof treatment and the acoustic and hygrothermal performances of the vegetal wools. After this study, it will be possible to expand the applications of vegetal wools for green building envelope.

Fig. 1 .

 1 (a) Hemp wool sample; (b) SEM observation of polymer and hemp fibers [23].

Fig. 2 .

 2 (a) Hemp wool sound absorption[START_REF] Piegay | Approche conjointe acoustique et thermique pour l'optimisation des laines végétales du bâtiment[END_REF]; (b) Hemp wool sorption curve[START_REF] Collet | Water vapor properties of two hemp wools manufactured with different treatments[END_REF]

Fig. 3 .

 3 Fig. 3. Methodology scheme

ACKNOWLEDGEMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie COFUND grant agreement No 101034248.