Juliette Legrand
email: ajuliette.legrand1@inrae.fr

Pierre Ailliot
email: cpierre.ailliot@univ-brest.fr

Philippe Naveau
email: philippe.naveau@lsce.ipsl.fr

Nicolas Raillard
email: nicolas.raillard@ifremer.fr

JOINT STOCHASTIC SIMULATION OF EXTREME COASTAL AND OFFSHORE SIGNIFICANT WAVE HEIGHTS

Keywords: Bivariate extremes, Multivariate generalised Pareto distribution, Simulation of extremes, Nonstationarity, Extended generalised Pareto distribution, Covariate effects, Significant wave heights. 1

The characterisation of future extreme wave events is crucial because of their multiple impacts, covering a broad range of topics such as coastal flood hazard, coastal erosion, reliability of offshore and coastal structures. The main goal of this paper is to propose and study a stochastic simulator that, given offshore conditions (peak direction Dp, peak period Tp and moderately high significant wave heights Hs), produces jointly offshore and coastal extreme Hs, a quantity measuring the wave severity and which represent a key feature in coastal risk analysis. For this purpose, we rely on bivariate Peaks over Threshold and a nonparametric simulation scheme of bivariate GPD is developed. From this joint simulator, a second generator is derived, allowing for conditional simulations of extreme Hs. Finally, to take into account nonstationarities, the extended generalised Pareto model is also adapted, letting the parameters vary with specific sea state parameters Tp and Dp. The performances of the two proposed generators are illustrated on simulated data and then applied to the simulation of new extreme oceanographic conditions close to the French Brittany coast using hindcast sea state data. Results show that the proposed algorithms successfully simulate future extreme Hs near the coast in a nonparametric way, jointly or conditionally on sea state parameters from a coarser model.

1. Introduction. French coastlines have been particularly affected by extreme maritime events in the past (Nicolae [START_REF] Lerma | Spatial variability of extreme wave height along the Atlantic and channel French coast[END_REF]. Co-occurrence of high tidal coefficients, atmospheric surge conditions and specific sea states can lead to extreme maritime events. These events are particularly crucial for assessing flooding risks and their consequences [START_REF] Genovese | Storm surge disaster risk management: the Xynthia case study in France[END_REF][START_REF] Bertin | Importance of wave age and resonance in storm surges: The case Xynthia, Bay of Biscay[END_REF]. According to the special IPCC report (see [START_REF] Collins | Extremes, Abrupt Changes and Managing Risk[END_REF]), extreme wave heights, which contribute to these extreme maritime events, have increased over the past few years. The recent IPCC report [START_REF] Seneviratne | Weather and Climate Extreme Events in a Changing Climate[END_REF] indicates, with high confidence, an increase in the occurrence and magnitude of such coastal events in the future. More specifically, [START_REF] Caires | 100-Year Return Value Estimates for Ocean Wind Speed and Significant Wave Height from the ERA-40 Data[END_REF] showed that the most extreme wave conditions were expected to occur in the North Atlantic, which includes the Bay of Biscay, our study area.

Sea surface elevation over a geographical area results from the superposition of waves generated by local winds and by remote swell (generated in distant regions). The characterisation of this complex surface is called a sea state and to describe it, various parameters are available. In this work, we will focus on three variables: the significant wave height denoted H s [m], the peak period, T p [s], and the peak direction, D p [°] (see e.g. [START_REF] Holthuijsen | Waves in Oceanic and Coastal Waters[END_REF] for more details). (left) Portion of IOWAGA hindcast database grid, the red dot corresponds to the "offshore" point (data extracted from the IOWAGA database) and the triangular green dot corresponds to the "coastal" point (data extracted from the HOMERE database). (right) Scatter plot of the coastal significant wave heights versus the offshore significant wave heights with the different thresholds considered. The dark dots belong to the region of the data used for simulation.

From a coastal risk point of view, a fundamental question is to determine how moderately high offshore significant wave heights can produce large coastal H s . Peak direction and peak period influence the relationship between coastal and offshore H s . In this context, our main goal is to propose and study a stochastic simulator that, given offshore conditions (T p , D p , H s moderately high), produces jointly offshore and coastal extreme significant wave heights. Note that in this study, we assume the availability of offshore T p and D p , but we could also consider an appropriate bivariate generator of (T p , D p ) (based, for example, on Heredia-Zavoni and [START_REF] Heredia-Zavoni | Modeling directional environmental contours using three dimensional vine copulas[END_REF]. Hereinafter, H c (resp. H o ) will denote the coastal (resp. offshore) significant wave heights. The left-hand side map of Figure 1 shows the two locations of interest.

From such a stochastic generator (first goal), many products can be derived. In particular our second objective is to propose a conditional simulation model (second goal). The framework for each step of this study is summarised in Table 1. In order to make the two simulation models as flexible as possible, nonparametric algorithms are derived using resampling techniques (or nonparametric bootstrap [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]). While this study focuses on simulation of extreme H s , and as illustrated by the numerical simulations in Section 4, the two nonparametric algorithms developed could be applied to a broad range of data.

In multivariate extreme value analysis, one is often interested in the joint behaviour of the variables as they become large. As illustrated by the right-hand side of Figure 1, which depicts a scatter plot between H o and H c , large values tend to occur simultaneously. For this Table 1: Summary of available data for each step of this study. A tick v (resp. a cross x)

indicates the availability (resp. non-availability) of the data.

Hc Ho Dp Tp

Inference v v v v First Goal: Joint simulation x x v v Second Goal: Conditional simulation x v v v
specific type of dependence, called asymptotic dependence [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], models from the class of multivariate Extreme Value Theory (EVT) can be used. To achieve the two objectives (first Goal and second Goal), we will therefore propose two simulation algorithms based on multivariate Peaks over Thresholds (Sec. 8.3. 1 Beirlant et al., 2004). Note that to assess whether data fall within the class of asymptotic dependence or not, summary statistics have been developed such as the dependence measures χ and χ [START_REF] Coles | Dependence Measures for Extreme Value Analyses[END_REF] and, for our specific dataset, these measures do not point toward the case of asymptotic independence between H o and H c (see Figure 11 in Appendix A).

For independent extremes, conditional models based on [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF] should be favoured to deal with our second goal, see for example [START_REF] Towe | Statistical downscaling for future extreme wave heights in the North Sea[END_REF]; [START_REF] Shooter | On spatial conditional extremes for ocean storm severity[END_REF]; [START_REF] Tendijck | Modeling the Extremes of Bivariate Mixture Distributions With Application to Oceanographic Data[END_REF]. In particular, if the distance between the two locations increases, the dependence is likely to decrease and other specific models should be used to address such issues (e.g. [START_REF] Shooter | On spatial conditional extremes for ocean storm severity[END_REF]Shooter et al., , 2021a)).

Before modelling the joint behaviour of large values, it is necessary to model margins [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. To visualise this task with respect to our data, the empirical histograms displayed in Figure 2 indicate that, given {H o > v o } (moderately high offshore significant wave heights), a traditional univariate extreme value approach based on fitting a generalised Pareto distribution (GPD) to the exceedances [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] is not appropriate. Indeed, the clear increase from zero to the mode in the left-hand histogram of Figure 2 cannot be reproduced by the probability density function (pdf) of the GPD that is a strictly decreasing function. To tackle this issue, we use the extended generalised Pareto distribution (EGPD) introduced by [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] which handles this type of setting (see also [START_REF] Papastathopoulos | Extended generalised Pareto models for tail estimation[END_REF]), more details are given in Section 3.

Furthermore, like many other environmental data, extreme H s are nonstationary with respect to covariates [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF] and marginal models need to take into account this nonstationarity (e.g. [START_REF] Ewans | The Effect of Directionality on Northern North Sea Extreme Wave Design Criteria[END_REF][START_REF] Méndez | Seasonality and duration in extreme value distributions of significant wave height[END_REF][START_REF] Casas-Prat | A physical-based statistical method for modeling ocean wave heights[END_REF]. To incorporate nonstationarities, [START_REF] Chavez-Demoulin | Generalized Additive Modelling of Sample Extremes[END_REF] proposed to let the parameters of an extreme value model vary as smooth functions of covariates. This has been intensively applied to oceanographic data (e.g. [START_REF] Feld | Estimation of Storm Peak and Intra-Storm Directional-Seasonal Design Conditions in the North Sea[END_REF]Fig 2: (left) Empirical histogram for the coastal significant wave height threshold exceedances and (right) similarly for the offshore significant wave heights threshold exceedances, illustrating that fitting a generalised Pareto distribution [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] is not suitable for the coastal data. The coastal marginal threshold is defined by

v c := min (H c ; H o > v o ).
Jonathan, Ewans and Randell, 2014;[START_REF] Ross | Efficient estimation of return value distributions from non-stationary marginal extreme value models using Bayesian inference[END_REF]. However, there are only a few papers dealing with nonstationary EGPD [START_REF] De | An Extreme Value Bayesian Lasso for the Conditional Left and Right Tails[END_REF][START_REF] Haruna | Performance-based comparison of regionalization methods to improve the at-site estimates of daily precipitation[END_REF]. In this study, and as illustrated in Figure 3, the marginal EGPD models parameters will be described as smooth functions of the covariates T p and D p . for different quantile levels q ∈ {0.25, 0.5, 0.75, 0.9, 0.99}, estimation is performed using smoothed quantile regression [START_REF] Koenker | Quantile smoothing splines[END_REF].

The key steps of our study are the following: (1) marginal regression modelling within the class of EGPD, (2) transformation of the data to common exponential margins, (3) modelling extremal dependence between the variables using multivariate generalised Pareto model (hereafter MGP models) [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF], (4) nonparametric simulation of bivariate extreme H s within the class of MGP distributions. In our modelling scheme, different steps are novelties.To our knowledge, little attention has been paid in the literature to the modelling of multivariate nonstationary extremes using EGPD and to the nonparametric simulation within the MGP class.

Our paper is organised as follows. In Section 2, the sea state data are presented and the marginal inference incorporating covariates in the EGPD modelling is described in Section 3. In Section 4, the nonparametric method to simulate MGP vectors is presented and some numerical experiments are shown. Two algorithms are outlined, one for bivariate simulations and a second one for conditional simulations. Both algorithms are applied in Section 5 to the sea state data. The R codes along with the data used in this study can be found in the Supplementary Material A [START_REF] Legrand | Supplement to[END_REF]).

2. Sea state data. Our study is carried out in the northern part of the Bay of Biscay in France. The specificity of this region is that it is exposed to the Atlantic Ocean and therefore subject to complex superpositions of wind generated waves and swell. The data are extracted from two different wave hindcasts provided by IFREMER and consist of simulations of sea states by a numerical model. First, the IOWAGA database [START_REF] Ardhuin | IOWAGA[END_REF] corresponds to sea states parameters that are generated by the wave model WAVEWATCH-III and forced by CFSR winds on a Global grid (0.5 • resolution grid in latitude and longitude). HOMERE is the second database [START_REF] Accensi | HOMERE[END_REF], also based on WAVEWATCH-III model and forced by IOWAGA on the wet boundaries, but on an unstructured grid covering only the English Channel and the Bay of Biscay, more refined close to the coast and with the inclusion of currents and water levels.

As mentioned in the introduction, we restrict our attention to two specific locations: an offshore grid point (47 • N, 3 • W) from the IOWAGA database and a coastal point from the HOMERE database, near the French coast (47 • 24N, 2 • 78W) corresponding to the SEM-REV sea test site [START_REF] Mouslim | Development of the French Wave Energy Test Site SEM-REV[END_REF] (see Figure 1). Among all the sea states parameters, the significant wave heights H s [m] from the two locations, the peak period T p [s] and the peak direction D p [ • ] only from the offshore location are used. Data are available at 3-hour intervals spanning from 1994 to 2016. More precisely, the original HOMERE database has a 1-hour resolution time step but the IOWAGA database is sampled every 3 hours, so to obtain data at the same time scale a sub-sampling of the HOMERE database every 3 hours is performed. Recall that H c (resp. H o ) correspond to the coastal (resp. offshore) significant wave heights. A scatter plot between H c and H o can be found in Figure 1, highlighting a strong dependence structure between the variables. Data are then split into two sets:

• Set 1 contains the first 70% of the data and is used for the inference of the marginal regression models and the preliminary steps for the simulation of H s (see Section 5); • Set 2 contains the remaining 30% of the data and is used for the simulation of extreme H s .

In order to propagate the uncertainties from the marginal modelling to the bivariate simulations, a bootstrap resampling is also adopted by repeating the combined marginal and bivariate analysis for 100 resamples of the Set 1. Note that Set 2 is not involved in the bootstrap procedure.

Marginal regression analysis.

3.1. Marginal regression. In this section, only the data from Set 1 are considered. A regression model for H c and H o is chosen. We pre-select the extremes by considering, within the Set 1, data such that H o > v o , i.e. belonging to the right rectangular region in Figure 1, where v o is defined as the 0.98 quantile of H o . A common choice when someone is interested in extreme values is to work with the class of the generalised Pareto distributions (GPD) (e.g. [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]). However, this type of model always raises questions on the choice of the threshold and the GPD approximation holds true only for the very high values. In our case we want to model all the data that are to the right of the vertical blue line in Figure 1, this means that values are not necessarily extremes. To overcome such problems, [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] proposed a new class of extreme value distributions, called extended generalised Pareto distributions (EGPD). The EGPD class is suitable for modelling the entire range of data, not only the most extreme values, and avoids the need for careful threshold selection. It also ensures that both lower and upper tails are in compliance with univariate EVT. In [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], four parametric models are proposed. We restrict ourselves to the first and simplest one (corresponding to the EGP3 model introduced by Papastathopoulos and Tawn ( 2013)), which appears to be flexible enough, and whose cumulative distribution function is of the form ( 1)

F (x) = 1 -1 + ξx σ -1/ξ κ .
The model has three parameters: scale σ > 0, shape ξ ∈ R and an additional parameter κ > 0 which controls the shape of the lower tail. In [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF], the authors developed the EGPD for non-negative shape parameter ξ. Indeed, the main applications were the modelling of daily rainfall (e.g. [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF][START_REF] Tencaliec | Flexible semiparametric generalized Pareto modeling of the entire range of rainfall amount[END_REF][START_REF] De | An Extreme Value Bayesian Lasso for the Conditional Left and Right Tails[END_REF][START_REF] Rivoire | A Comparison of Moderate and Extreme ERA-5 Daily Precipitation With Two Observational Data Sets[END_REF], which are heavy-tailed (ξ > 0). In our case, as reported in [START_REF] Ewans | Statistical modelling of extreme ocean environments for marine design: A review[END_REF], extreme H s data are generally described by upperbounded tail distributions. Still, the case ξ < 0 can be handled by model [START_REF] Casas-Prat | A physical-based statistical method for modeling ocean wave heights[END_REF]. As mentioned in the Introduction, H s data are nonstationary (see Figure 3), see also [START_REF] Jonathan | Non-stationary conditional extremes of northern North Sea storm characteristics[END_REF]; [START_REF] Feld | Estimation of Storm Peak and Intra-Storm Directional-Seasonal Design Conditions in the North Sea[END_REF][START_REF] De | Non-stationary extreme value analysis of sea states based on linear trends. Analysis of annual maxima series of significant wave height and peak period in the Mediterranean Sea[END_REF]. Therefore, we regress H s on the peak direction D p and the peak period T p . We choose here to put the dependency on the scale parameter. The regression marginal models can then be written as follows

(2)            P(H o -v o ≤ x|H o > v o , T p , D p ) = 1 -1 + ξ o x σ o (T p , D p ) -1/ξo κo , P(H c -v c ≤ x|H o > v o , T p , D p ) = 1 -1 + ξ c x σ c (T p , D p ) -1/ξc κc . The coastal marginal threshold v c introduced in (2) is defined by v c := min(H c ; H o > v o ) so that the minimum of H c -v c is equal to zero.
This definition of v c allows to keep a large amount of data for the simulations (see Section 5) and the use of the EGPD ensures that extreme values of H c are still well modelled.

Note that in classical EVT, incorporating a varying threshold in a extreme value model often leads to improved marginal modelling (see for example [START_REF] Northrop | Threshold modelling of spatially dependent nonstationary extremes with application to hurricane-induced wave heights[END_REF] for a spatially varying threshold to model extreme H s ). However, in Equation ( 2), the thresholds v c and v o should not be confused with the usual EVT thresholds. For example, going back to the left panel of Figure 2, the histogram clearly highlights a mode around 1.5. In this case, a classical EVT approach would have probably been based on varying thresholds greater than 1.5. The parameter κ of the EGPD allows us to capture the lower tail behaviour. This additional parameter appears to bring the necessary flexibility to fit well this particular dataset, see Figure 12, without the need of a varying EVT threshold. If needed, one could integrate covariate in κ (Le Carrer, 2022).

The regression marginal models (2) are estimated using the R package gamlss [START_REF] Stasinopoulos | Instructions on how to use the gamlss package in R Second Edition[END_REF] with the EGPD family (Le Carrer, 2022). The inference is performed using maximum penalised likelihood estimation (note that the model fitting is achieved with the CG algorithm [START_REF] Cole | Smoothing reference centile curves: The LMS method and penalized likelihood[END_REF]). In our model (2), we assume that the parameters σ o and σ c vary smoothly with T p and D p . This is achieved using tensor product of cubic regression splines. The parameter estimates for the regression marginal models are reported in Table 2. Both asymptotic 95% confidence intervals (CI), derived from the asymptotic variance-covariance matrix of the fitted models, and bootstrap 95% CI for each parameter are given in brackets. Both procedures produce similar results.

Both estimated shape parameters are negative but close to zero, which suggests lighttailed or bounded distributions. This is in accordance with previous studies and the physical behaviour of wave heights in shallow waters [START_REF] Castillo | Engineering Analysis of Extreme Value Data: Selection of Models[END_REF]Vanem and Fazeres-Ferradosa, 2022). Regarding the lower tail parameter κ, looking at Equation ( 1) one can check that less mass is put at zero for large values of κ. Thus looking at the empirical histograms in Figure 2, the estimates of κ seem reasonable. Goodness-of-fit plots of the model (2) are shown in Appendix B.

Covariate effects.

The covariate effects on the scale parameters can be seen on the top panels in Figure 4. Regarding the effect of D p , we find that the biggest storms come from the WSW and with strong similarities between the coastal and the offshore model. There are more differences between the two locations if we look at the effect of T p , which might be due to a greater loss of energy during the propagation of waves with large period. This is in particular the case for swells generated well offshore, corresponding to a large peak period, and coming from the NW, which are filtered at the coast due to the islands. To visualise in an alternative way the effect of the covariates on the scale parameters, we also choose to consider the theoretical expectation of the fitted EGPD models. As can be seen from Equation (3) (and similarly for H o ), the theoretical expectation of model ( 2) is directly proportional to the scale parameter (see [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF]):

(3)

E(H c |H o > v o , T p , D p ) = σ c (T p , D p ) 1 ξ c [κ c B(κ c , 1 -ξ c ) -1] ,
where B denotes the Beta function defined by The bottom left panel in Figure 4 represents the empirical ratio H c /H o given T p and D p values. Local polynomial interpolation (LOESS, [START_REF] Cleveland | Locally Weighted Regression: An Approach to Regression Analysis by Local Fitting[END_REF]) is performed between the observed data points to get values on a regular grid of D p and T p , ranging from 230 to 300 degrees for the peak direction and from 10 to 20 seconds for the peak period. This first panel is then compared to the estimated ratio of the two conditional expectations (bottom right panel in Figure 4), which, from Equation (3), is proportional to

B(a, b) = 1 0 t a-1 (1 -t) b-1 dt.
σc (T p , D p ) /σ o (T p , D p ) .
This ratio can give us an idea of the propagation of the wave energy from the offshore to the coast as a function of the covariates. These results can be physically interpreted: the loss of wave energy between the offshore and the coast is lower for small periods but also for waves coming from the SW rather than the NW due to the bathymetry (see the map on Figure 1).

Using the estimated σc (T p , D p ) and σo (T p , D p ), the H s data are then transformed to common exponential scale using the probability integral transform ( 4)

H E o := -log 1 -Fo [(H o -v o )/σ o (T p , D p )] , H E c := -log 1 -Fc [(H c -v c )/σ c (T p , D p )] ,
where Fo (resp. Fc ) is the fitted EGP D( ξo , κo , 1) cdf's (resp. EGP D( ξc , κc , 1)) from Equation (2).

4. Multivariate Pareto model. In this section, the threshold exceedances of H s transformed to common exponential margins are modelled. This vector is denoted

H E := (H E c , H E o ) in the following.
For that, we adapt the definition of [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] of bivariate threshold exceedances given as ( 5)

H E -u|H E ≰ u where u := (u c , u o ) ∈ R 2 + and H E ≰ u means that H E c > u c and/or H E o > u o
, that is to say we are extreme in at least one of the two components. Then multivariate EVT theory states that (5) can be well approximated by a multivariate generalised Pareto (MGP) distribution [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF]. Note that there are different equivalent definitions for multivariate threshold exceedances (see Section 8.3 of [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]).

Rootzén, Segers and Wadsworth (2018) derived a stochastic representation of standard MGP vectors considering that a bivariate random vector Z follows a standard MGP distribution if, and only if,

(6) Z = E + T -max(T ),
with T a random vector and E a unit exponential random variable independent of T . In the above equation, the addition should be interpreted component-wise and the scalars E and max(T ) are repeated twice if Z is a bivariate random vector. Note that the standard MGP distribution is supported by the set L := {x ∈ R d ; x ≰ 0}.

In our study, Equation ( 6) is adapted to Z = (Z 1 , Z 2 ) defined as ( 7)

Z 1 := H E o -u o |H E o > u o or H E c > u c , Z 2 := H E c -u c |H E o > u o or H E c > u c . Algorithm 1 Nonparametric bootstrap MGP simulation 1: input A sample (Z 1,i , Z 2,i ) 1≤i≤n from a MGP distribution 2: output A simulated sample (Z (m) 1,k , Z (m) 2,k ) 1≤k≤m , potentially with m ̸ = n 3: procedure 4: Define ∆ i := Z 1,i -Z 2,i for 1 ≤ i ≤ n 5: Generate m realisations E (m) k ∼ Exp(1), independently of (∆ i ) 1≤i≤n , for 1 ≤ k ≤ m 6: Bootstrap m realisations ∆ (m) k , 1 ≤ k ≤ m, from (∆ 1 , . . . , ∆n) 7: end procedure 8: return Z (m) 1,k := E (m) k + ∆ (m) k 1 ∆ (m) k <0 and Z (m) 2,k := E (m) k -∆ (m) k 1 ∆ (m) k >0
, for 1 ≤ k ≤ m 4.1. Simulation of bivariate standard generalised Pareto distributed vectors. [START_REF] Kiriliouk | Peaks Over Thresholds Modeling With Multivariate Generalized Pareto Distributions[END_REF] established several parametric MGP models by setting explicit densities for T in a multivariate setting. In the following, we consider only vectors of dimension 2, i.e. Z = (Z 1 , Z 2 ) and T = (T 1 , T 2 ). To bypass the choice of the underlying distribution for T , we start from the following rewriting of Equation ( 6), ( 8)

Z 1 = E + ∆1 ∆<0 , Z 2 = E -∆1 ∆≥0 ,
where ∆ := Z 1 -Z 2 = T 1 -T 2 and 1 A denotes the indicator function, equals to 1 if A is true and 0 otherwise. Equation ( 8) is the basis for our simulation algorithms. From this equation, we see that we need to simulate values of ∆ and E independently, instead of (T 1 , T 2 ). Generating independent, and identically distributed, unit exponentials is trivial, so the main difficulty is to simulate ∆. This can be achieved by bootstrapping (see [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]). Our approach is then described in Algorithm 1 and a theoretical proof can be found in Appendix C.

Numerical experiments.

In the following, we simulate MGP vectors Z = (Z 1 , Z 2 ) from the representation (6) with different parametric models on (T 1 , T 2 ) and we compare with our simulation algorithm. The different experiments are reported in Table 3 and some graphical results are shown in Figure 5 which displays for each model a scatter plot of the data, the measure of extremal dependence χ(u) for increasing values of u, and the marginal quantile-quantile plots. We use the measure χ(u) which gives a measure of asymptotic dependence between two variables X and Y (for more details see e.g. Coles, [START_REF] Coles | Dependence Measures for Extreme Value Analyses[END_REF]) and which is defined by ( 9)

χ(u) := P Y > F -1 Y (u) | X > F -1 X (u)
, u ∈ (0, 1). In Table 3, (a) and (b) are two bivariate Gaussian models with same correlation coefficient ρ < 1 but with µ 1 ̸ = µ 2 for (b), leading to asymmetry. Model (c) corresponds to the Type I bivariate logistic distribution proposed by [START_REF] Gumbel | Bivariate Logistic Distributions[END_REF]. For model (d) we consider two independent Gumbel distributed variables with different scale parameters. And lastly, (e) corresponds to a bivariate exponential distribution as defined in [START_REF] Marshall | A Multivariate Exponential Distribution[END_REF].

As seen in Figure 5, Algorithm 1 successfully simulates draws from the parametric simulations in terms of the marginals Z 1 and Z 2 , but also recovers well the dependence structure when looking at the measure of dependence χ(u). Table 3: Overview of the different experiments carried-out given the joint distribution of T .

For each, we give the joint distribution F (x 1 , x 2 ) when it can be written easily or the survival function S(x 1 , x 2 ). In the third column, we give the different parameters values used in the numerical experiments.

Bivariate model

Joint distribution of T Parameters (a) Gaussian symmetric N (µ 1 , µ 2 ), 1 ρ ρ 1 µ 1 = 0 µ 2 = 0 ρ = 0.4 (b) Gaussian asymmetric N (µ 1 , µ 2 ), 1 ρ ρ 1 µ 1 = 0 µ 2 = 2 ρ = 0.4 (c) Logistic F (x 1 , x 2 ) = 1 + e -x 1 /σ 1 + e -x 2 /σ 2 -1 σ 1 = 1 σ 2 = 5 (d) Gumbel F (x 1 , x 2 ) = exp [-exp{-x 1 /σ 1 }] exp [-exp{-x 2 /σ 2 }] σ 1 = 1 σ 2 = 4 (e) Exponential S(x 1 , x 2 ) = exp {-λ 1 x 1 -λ 2 x 2 -λ 3 max(x 1 , x 2 )} λ 1 = 2 λ 2 = 10 λ 3 = 1
4.3. Conditional simulation within the MGP class. From an application perspective, we also want to be able to simulate conditionally on one of the two variables. In this section we describe the conditional simulation algorithm for the MGP model.

From Equation ( 6) we can derive (10)

Z 2 = Z 1 + T 2 -T 1 = Z 1 -∆.
From Equation (10), we can design a simulation strategy but caution is required because ∆ and Z 1 are not necessarily independent. But for some values of Z 1 , this will be the case.

To see this, one can compute the conditional distribution of ∆ given Z 1 = z 1 starting from the joint distribution function of (Z 1 , Z 2 ) which is given by

f (Z1,Z2) (z 1 , z 2 ) = e -max(z1,z2) f ∆ (z 1 -z 2 )1 (z1,z2)∈L , for z 1 , z 2 ∈ R,
where f ∆ denotes the distribution function of ∆.

First case: If z 1 > 0. In this case, noting that if z 1 > 0 then 1 (z1,z2)∈L = 1, the marginal distribution of Z 1 is given as follows

f Z1 (z 1 ) = e -z1 K, where K := 0 -∞ e u f ∆ (u)du + ∞ 0 f ∆ (u)du and does not depend on z 1 . Therefore, the conditional distribution of Z 2 given Z 1 = z 1 when z 1 > 0 is given by f Z2|Z1 (z 2 | z 1 ) = 1 K f ∆ (z 1 -z 2 )1 z1≥z2 + e z1-z2 f ∆ (z 1 -z 2 )1 z1<z2
From this, the conditional distribution of ∆ given Z 1 = z 1 > 0 can then be deduced:

(11) f ∆|Z1 (δ | z 1 ) = 1 K f ∆ (δ)1 δ≥0 + e δ f ∆ (δ)1 δ<0 .
This shows that, conditionally on Z 1 > 0, ∆ does not depend on Z 1 .

Second case: If z 1 < 0. Then, noting that 1 (z1,z2)∈L = 1 z2>0 if z 1 < 0, we get

f Z1 (z 1 ) = e -z1 K(z 1 ), Algorithm 2 Nonparametric conditional MGP simulation 1: input A sample (∆ i ) 1≤i≤n ; a realisation z 1 of Z 1 2: output A simulated sample (Z (m) 2,k ) 1≤k≤m conditionally on Z 1 = z 1 , potentially with m ̸ = n 3: procedure 4: if z 1 > 0 then 5: Define ∆ |Z + 1 the subset of (∆ i ) 1≤i≤n such that Z 1 > 0 6: Bootstrap m realisations ∆ (m) k , 1 ≤ k ≤ m, from ∆ |Z + 1 independently of Z 1 7: else 8: for 1 ≤ k ≤ m do 9: Sample one realisation ∆ (m) k from (∆ i ) 1≤i≤n independently of Z 1 10: Generate a random number u ∈ [0, 1] 11: while u > exp(∆ (m) k )1 ∆ (m) k <z 1 do 12:
Repeat steps 9 and 10 13:

end for 14: end procedure

15: return Z (m) 2,k := z 1 -∆ (m) k for 1 ≤ k ≤ m where K(z 1 ) := z1 -∞ e u f ∆ (u)du. And we can derive the conditional distribution of Z 2 given Z 1 = z 1 < 0 as follows f Z2|Z1 (z 2 | z 1 ) = 1 K(z 1 ) e z1-z2 f ∆ (z 1 -z 2 )1 z2>0 .
The conditional distribution of ∆ given Z 1 = z 1 < 0 is then given by ( 12)

f ∆|Z1 (δ | z 1 ) = 1 K(z 1 )
e δ f ∆ (δ)1 δ<z1 .

From Equations ( 11) and ( 12) we derive the conditional simulation algorithm described in Algorithm 2, where the simulation procedure is split into two cases:

1. If z 1 > 0, we can sample values of ∆ independently of Z 1 , 2. otherwise, if z 1 < 0, we use a rejection sampling approach to approximate the targeted conditional density in Equation ( 12).

Similarly, we could also derive a simulation scheme of Z 1 given Z 2 = z 2 .

4.4. Numerical experiment continued. As for the bivariate simulations, we can illustrate Algorithm 2 with numerical experiments. We choose here to show the results only for Model (a) (Symmetric Gaussian) since for this specific model we have an explicit form for the theoretical distribution of ∆. The results are presented in Figure 6 where we simulated the conditional distribution of Z 2 for eight different conditioning values. The sampled and theoretical conditional distributions appear to be in close conformity.

5. Application to extreme significant wave height. The methodology presented in Section 4 is applied to the joint and the conditional simulations of extreme significant wave heights. For that, the sample of bivariate threshold exceedances (Z 1 , Z 2 ) defined in Equation ( 7) is used as input data for Algorithm 1 or Algorithm 2. The thresholds u o and u c in (7) are defined as the 0.8 quantile of H E o , or equivalently of H E c . This quantile has been selected following the stability plot approach proposed in [START_REF] Kiriliouk | Peaks Over Thresholds Modeling With Multivariate Generalized Pareto Distributions[END_REF].

Recall that in both cases, simulations are performed on the exponential scale. A final step of back transformation is then necessary to get simulations of H s on the original scale. This final step corresponds to part 3 (resp. 5) in the following procedure for the joint (resp. conditional) simulation of H s . For the sake of clarity we now divide the joint and the conditional simulation scheme of H s in two separate sections. Note that in the following two sections, simulations are performed given the marginal parameters estimates. The quantification of the uncertainty propagation from the marginal inference step to the bivariate estimation is postponed to Section 5.3. 5.1. Joint simulation of significant wave heights. The joint simulation scheme for extreme H s is described hereafter. In the following we fix the pair value (t p , d p ) ∈ R 2 which may be taken from Set 2.

1. Compute σo (t p , d p ) and σc (t p , d p ) from the marginal EGPD models fitted on Set 1 (see Section 3). 2. Simulate m pairs of (z 1 , z 2 ) applying Algorithm 1 with input data (Z 1 , Z 2 ) as defined in (7). We therefore obtain m simulated pairs ((z 1,1 , z 2,1 ), . . . , (z 1,m , z 2,m )) for a fixed value (t p , d p ).

Transform the simulated values to the original scale

h o,i := σo (t p , d p ) F -1 o (1 -e -(z1,i+uo) ) + v o ∈ R m , h c,i := σc (t p , d p ) F -1 c (1 -e -(z2,i+uc) ) + v c ∈ R m ,
where

F -1 o (resp. F -1 c ) is the inverse cdf of the EGP D( ξo , κo , 1) (resp. EGP D( ξc , κc , 1)) estimated in Section 3.
This procedure is then applied to four selected pairs (t p , d p ) from the Set 2 corresponding to the four largest H c of this dataset. Figure 7 depicts simulated pairs of offshore and coastal H s with simulation sample size m = 1000.

Figure 7 shows that for these specific conditions, Algorithm 1 successfully generates extreme H c and H o . Note that since the four points considered are the largest observations, they are expected to be among the extremes of the simulated distributions. 4 gives the empirical estimate of the joint survival probability defined as 1/m m i=1 1{h o,i > h o , h c,i > h c }, applying the above procedure with simulation sample size m = 1.10 6 . This estimated probability quantifies the observations made with Figure 7: for the most extreme events, the associated joint probabilities are expected to be lower. 5.2. Conditional simulation of coastal significant wave heights. The conditional simulation scheme for extreme coastal H s is described below. Note that the procedure is symmetrical for simulating offshore H s . In the following we fix the triplet value (h o , t p , d p ) ∈ R 3 which may be taken from Set 2. Table 4: Empirical joint survival probability of exceeding the observed extreme significant wave heights h c and h o for the eight largest coastal significant wave heights of Set 2. Only events from different storms are given (i.e. events separated by more than 3 hours). Estimation is performed using m = 10 6 simulated pairs (h o,i , h c,i ) 1≤i≤m for each largest event with Algorithm 1.

Date-Time Hc

[m] Ho [m] Tp [s] Dp [°] Joint
h E o = -log 1 -Fo [(h o -v o )/σ o (t p , d p )]
where Fo is the EGP ( ξo , κo , 1) cdf's. 3. Set z 1 := h E o -u o , where u o is the threshold on the offshore H s on the exponential scale. 4. Simulate z 2 applying Algorithm 2, only in the case z 1 > 0. Here ∆ |Z + 1 is defined from Set 1 through ∆ |Z + 1 := (Z 1 -Z 2 )1{Z 1 > 0}, with Z 1 , Z 2 as defined in ( 7), and bootstrapped m times. We therefore obtain m simulations of z 2 = (z 2,1 , . . . , z 2,m ) for a fixed triplet (h o , t p , d p ), given z 1 > 0.

Transform the predicted values to the original scale

h c,i := σc (t p , d p ) F -1 c (1 -e -(z2,i+uc) ) + v c where F -1 c
is the inverse cdf of the EGP D( ξc , κc , 1) (see Section 3).

Note that in

Step 4 above, the simulation is restricted to the case when z 1 > 0 for convenience, since our focus is on the simulation of extreme H s .

The pseudo-algorithm described above is applied with all the triplet values (h o , t p , d p ) from Set 2, with simulation sample size m = 1000. These conditional simulations of coastal significant wave heights are then compared to the true values of H c from Set 2. The overall coverage probability (i.e. the number of times the actual value of coastal H s is within the 95% range of the predicted distribution) is equal to 95% and the simulations are shown in Figure 8. The simulations and the true H c values (red dots) are most of the time in good agreement. Since no declustering approach has been adopted, consecutive observations, which belong to the same storm event, are kept. They are depicted with identical colour, highlighting a temporal dependence structure between each storm cluster.

The effect of the covariates T p and D p in the conditional simulations is depicted in Figure 9, showing that the simulation model is able to simulate both the most intense and the more moderate H c . This plot also highlights for which sea state conditions the simulations are far from the observed values. It appears that the two predictions such that the observed H c value does not fall within the 95% simulation range correspond to small H c and H o .

Uncertainty quantification.

As mentioned in Section 2, in addition to these simulations, we also combined the marginal and the dependence analysis in a bootstrap scheme in order to accommodate for marginal uncertainties. Both Algorithm 1 and Algorithm 2 have been applied within the bootstrap procedure, leading to 100 bootstrap sample of joint and conditional simulations of H s . Results of each iteration can be found in the Supplementary Material B [START_REF] Legrand | Supplement to[END_REF]), showing the reproduction of Figures 7 and8 for each bootstrap sample. Overall, it seems that the performances are reasonably good within each resampling iteration. Also, regarding the conditional simulation algorithm, the mean coverage probability over the 100 bootstraps is equal to 93%.

6. Discussion and conclusions. Simulation of extreme events in a multivariate setting is of great interest to capture not only the statistical behaviour of the extremes, but also the dependence between large values of complex processes. Based on the multivariate EVT, this work presents two nonparametric simulation algorithms of bivariate generalised Pareto distributed variables, without assuming any specific parametric shape for the MGP model. Thanks to Algorithm 1, one can simulate joint extremes. As for Algorithm 2, it allows the simulation of conditional extremes. Both methods have been validated with numerical simulations.

We would like to point out that in the context of bivariate extremes, other simulation algorithms have been developed. For example [START_REF] Marcon | A semi-parametric stochastic generator for bivariate extreme events[END_REF] proposed a simulation method with a semi-parametric structure for the extremal dependence function, but it was not based on a MGP model and did not cover nonstationarities. [START_REF] Michel | Simulation and estimation in multivariate generalized Pareto models[END_REF] derived a nonparametric simulation framework of bivariate generalised Pareto variables using a different representation of MGP vectors than the one used in this paper.

For application purposes and as a by-product, a nonstationary marginal modelling with the EGPD was also developed, adding covariate effects on the scale parameter of the EGPD using smoothing splines.

We applied this work to the simulation of extreme significant wave heights near the Brittany coast given specific offshore sea state conditions (T p , D p ) with compelling results. In both joint and conditional settings, thanks to the presented algorithms, we are able to simulate realistic extreme H s events. Regarding the marginal models, other studies have considered extreme value models for H s with both varying scale and shape parameters (e.g. [START_REF] Ewans | The effect of directionality on extreme wave design criteria[END_REF][START_REF] Feld | Estimation of Storm Peak and Intra-Storm Directional-Seasonal Design Conditions in the North Sea[END_REF]. Incorporating effects of T p and D p on both parameters could certainly improve the models, but at the cost of computational limitations. Moreover, one could argue that for our specific dataset, the sea states could be considered as homogeneous since D p has a constrained domain between [234 • , 300 • ].

Note that in possible extensions of this work to climate projections, it is assumed that data will not be available at the coastal location but only on a coarse grid, similar to the IOWAGA Global hindcast. This argument, which is illustrated in Table 1, motivates the first pre-selection of the H s data through {H o > v o } for the marginal regression analysis (see Section 3).

Extensions to the multivariate case will be the subject of future works. Considering more than two locations raises different modelling issues. It would also be interesting to apply this methodology to other locations in order to ensure the proper generalisation of the methodology.

Finally, as already mentioned, a classical approach regarding extreme significant wave height simulations is the conditional extreme model of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF]. For example, [START_REF] Shooter | Spatial dependence of extreme seas in the North East Atlantic from satellite altimeter measurements[END_REF][START_REF] Shooter | Multivariate spatial conditional extremes for extreme ocean environments[END_REF] applied this model in a spatial setting to characterise the behaviour of the extremal dependence structure between different metocean variables. One strength of our bivariate model (hereafter referred to as MGPD approach) is that a conditional sampling scheme can be easily produced from simulations of joint extremes. So, we can compare the MGPD approach to the traditional conditional model of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF] for our data at hand. We consider hereinafter the conditional generation of extreme H s on the standardised scales. In other words, we generated standardised coastal significant wave heights on the exponential scale with Algorithm 2 on the one hand and with the conditional extreme value model of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF] on the other hand. We simulated n sim = 1000 samples of same length as the observation exceedances h E o > u o . Figure 10 depicts one of these n sim simulations of coastal H s versus observed offshore H s in the exponential space. One can see that for this specific simulated sample, the conditional extreme approach tends to simulate less extreme values of H c than the MGPD approach. A first comparison has been made computing the mean absolute error between the observed data and the simulations for both models. Results are depicted in the right panel of Figure 10, showing slightly better results for the MGPD model. 9)) and χ(u) between the offshore and the coastal significant wave heights. Recall that the dependence measure χ(u) (e.g. [START_REF] Coles | Dependence Measures for Extreme Value Analyses[END_REF] between two variables X and Y is define as

χ(u) = 2 log (1 -u) log P (F X (X) > u, F Y (Y ) > u)
, u ∈ (0, 1).

APPENDIX B: MARGINAL REGRESSION MODELLING

We show here the goodness of fit for the marginal regression models defined in Equation (2). As our models depend on some covariates, the diagnostic plots presented here are built for the standardised H s exceedances which are defined as (H c -v c ) /σ c (T p , D p ) (and similarly for H o ). From Figure 12, one can see that the fits seem to be fairly good, a slight discrepancy 2 ) ≤ (z 1 , z 2 ) = 1 -E min 1, e -min(z1-∆ (m) 1 ∆ (m) <0 ,z2+∆ (m) 1 ∆ (m) >0 ) = 1 -E min 1, e -min(z1-∆ (m) ,z2)-max(∆ (m) ,0) , for any (z 1 , z 2 ) ∈ {x ∈ R 2 ; x ≰ 0}.

Then, one can show that the function x → min (1, x), defined for x ≥ 0, is Lipschitz and bounded by 1. And applying the Portmanteau theorem, we have, letting min(n, m) → ∞,

P (Z (m) 1 , Z (m)
2 ) ≤ (z 1 , z 2 ) | ∆ 1 , . . . , ∆ n → 1 -E min 1, e -min(z1-∆,z2)-max(∆,0) = 1 -E min 1, e max(T1-z1,T2-z2)-max(T1,T2) . This is the cumulative distribution function of the MGP vector (Z 1 , Z 2 ) as defined in Rootzén, Segers and Wadsworth (2018) (Prop. 8).

The assumption in Lemma 1 is linked to bootstrap asymptotic theory (e.g. [START_REF] Bickel | Some Asymptotic Theory for the Bootstrap[END_REF]).

Fig 1 :

 1 Fig 1:(left) Portion of IOWAGA hindcast database grid, the red dot corresponds to the "offshore" point (data extracted from the IOWAGA database) and the triangular green dot corresponds to the "coastal" point (data extracted from the HOMERE database). (right) Scatter plot of the coastal significant wave heights versus the offshore significant wave heights with the different thresholds considered. The dark dots belong to the region of the data used for simulation.

Fig 3 :

 3 Fig 3: Dependence of H s with respect to T p and D p (left) Offshore significant wave heights H o given peak period T p . (right) Estimated quantiles of H o given the peak direction D pfor different quantile levels q ∈ {0.25, 0.5, 0.75, 0.9, 0.99}, estimation is performed using smoothed quantile regression[START_REF] Koenker | Quantile smoothing splines[END_REF].

Fig 4 :

 4 Fig 4: (top left) Estimated offshore scale parameter conditionally to the offshore peak direction D p and the offshore peak period T p . (top right) Identical to the top left panel but for the coastal scale parameter. (bottom left) Interpolated ratio of empirical extreme coastal significant wave heights H c and extreme offshore significant wave heights H o . The interpolated surface between the data points is performed using local polynomial interpolation of degree 2. (bottom right) Ratio of the predicted conditional expectations of extreme coastal significant wave heights H c and extreme offshore significant wave heights H o , conditionally to the offshore peak direction D p and the offshore peak period T p . On the four plots, observed data points are superimposed.

Fig 5 :

 5 Fig 5: Each panel line corresponds to one of the parametric model (a) to (e), and shows from left to right: (1) Scatter plot of simulated data with the parametric model with sample size n = 10000 (black dots) and sampled data from one simulation using Algorithm 1 with sample size m = 10000 (blue diamond-shaped dots); (2) Empirical estimates of the measure of asymptotic dependence χ(u) for the simulated data with the parametric model (black line) and for the sampled data from Algorithm 1 (red line), with associated 95% pointwise confidence intervals based on 1000 bootstrap replications; (3) and (4) Quantile-quantile plots for Z 1 and Z 2 with associated 95% pointwise confidence intervals based on 1000 bootstrap replications.

Fig 6 :

 6 Fig 6: Sampled conditional distribution of Z 2 given Z 1 = z 1 using Algorithm 2 for the asymmetric Gaussian model with sample size n = 10000. Eight experiments are presented for different quantiles of Z 1 whose values are reported in each panel title. The sample size for each simulation is m = 10000. The theoretical conditional density is superimposed in red.

Fig 7 :

 7 Fig 7: Sampled values of coastal versus offshore significant wave heights from Algorithm 1, m = 1000 pairs of points are simulated. Each panel corresponds to a fixed value of (t p , d p ) corresponding to the four largest H c from Set 2 in decreasing order (recall that Set 2 contains the last 30% of the original dataset that have not yet been used). Next to each scatter plot, marginal distribution of simulated H o and H c are depicted with boxplots. On both the scatter plots and the boxplots, the red diamond-shaped points represent the true values of the coastal and offshore significant wave heights for the four pairs.

Fig 8 :

 8 Fig 8: Boxplot of predicted H c conditionally on (H o , T p , D p ) using Algorithm 2. The simulation sample size is equal to m = 1000 for each observation. Red dots represents the observed H c values from Set 2. The alternating colours depict different storms: consecutive boxes with same colour correspond to observations that belong to the same storm (i.e. separated by less than 3 hours).

Fig 9 :

 9 Fig 9: Scatter plot from Set 2 with the peak direction D p on the y-axis and the peak period T p on the x-axis. A dot's colour corresponds to the value of the coastal significant wave heights H c . The size of the dots corresponds to the value of the offshore significant wave heights H o . The shape indicates if the observed H c falls within the 95% range of the predicted distribution from the conditional model where the simulation sample size is equal to m = 1000 for each observation.

Fig 10 :

 10 Fig 10: Scatter plots between coastal H s and offshore H s in the standardised space (exponential margins). The coloured triangle shaped points correspond to one simulation run with the conditional extreme model of Heffernan and Tawn (2004) (left) and with our MGPD approach (middle). The vertical dashed line corresponds to u o . (right) Mean absolute error over the n sim = 1000 simulations for both models: Heffernan and Tawn (2004) on the left -MGPD approach on the right.

Fig 11 :

 11 Fig 11: Empirical estimates of the measure of asymptotic dependence χ(u) (left) and χ(u) (right) between H o and H c , with 95% pointwise confidence intervals.

Fig 12 :

 12 Fig 12: (left) Empirical histogram (grey) of the standardised extreme significant wave heights exceedances (a) at the coast and (b) offshore. The fitted EGPD density is superimposed.(right) The corresponding quantile-quantile plots with associated 95% pointwise confidence intervals computed using parametric bootstrap.

Table 2 :

 2 Estimated parameters for the regression marginal models. Point estimates and 95% asymptotic and bootstrap confidence intervals are given in brackets.

	Parameter Estimate	Asymptotic CI	Bootstrap CI
	ξc	-0.11	[-0.15, -0.07] [-0.21, -0.09]
	κc	4.11	[3.57, 4.64]	[3.02, 4.84]
	ξo	-0.10	[-0.16, -0.04] [-0.17, -0.04]
	κo	1.16	[1.05, 1.26]	[1.05, 1.28]

  Table 4 describes the eight largest h c of Set 2, giving the time event and the corresponding h o , t p and d p values. The last column of Table

  Compute σo (t p , d p ) and σc (t p , d p ) from the marginal EGPD models fitted on Set 1 (see Section 3). 2. Transform h o to the standardised space using the probability integral transform:

	1.
						probability
	1 2014-02-05 12:00:00 GMT	7.83	10.1	17.86	257	0.002
	2 2011-12-16 03:00:00 GMT	7.74	10.4	13.89	262	0.045
	3 2010-02-24 15:00:00 GMT	7.62	10.0	15.15	259	0.057
	4 2016-01-02 06:00:00 GMT	7.18	8.4	14.71	256	0.277
	5 2013-12-24 06:00:00 GMT	7.08	8.4	14.71	253	0.388
	6 2014-02-14 21:00:00 GMT	7.06	8.3	14.08	248	0.691
	7 2015-01-15 09:00:00 GMT	6.45	7.9	13.70	260	0.489
	8 2016-03-28 03:00:00 GMT	6.30	7.6	14.08	256	0.762

VANEM, E. and FAZERES-FERRADOSA, T. (2022). A truncated, translated Weibull distribution for shallow water sea states. Coastal Engineering 172 104077. https://doi.org/10.1016/j.coastaleng.2021.104077

Acknowledgments. The authors would like to thank Noémie Le Carrer for helpful discussion on the implementation of the EGPD family within the gamlss package.

The authors are indebted to the Laboratoire d'Océanographie Physique et Spatiale, IFRE-MER for providing the data used in this study. Funding. The third author was supported in part by the DAMOCLES-COST-ACTION on compound events, the French national program (FRAISE-LEFE/INSU and 80 PRIME CNRS-INSU), and the European H2020 XAIDA (Grant agreement ID: 101003469. The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-20-CE40-0025-01 (T-REX project), the ANR-Melody (ANR-19-CE46-0011).

SUPPLEMENTARY MATERIAL

Supplement A: R code and data This supplement contains the R code and the data used in this study.

Supplement B: Bootstrap simulations

This supplement file contains two movies showing the results of each n boot = 100 bootstrap for both algorithms iteration by reproducing Figure 7 and Figure 8 of the main paper.


Brunovsky decomposition for dynamic interval localization

Simon Rohou 1 , Luc Jaulin 1 1 ENSTA Bretagne, Lab-STICC, UMR CNRS 6285, Brest, France Abstract-This paper proposes a new set-membership method for estimating the trajectories of dynamical systems, when the states are completely unknown and only non-linear observations are available. The first part of the proposed method is symbolic and follows the decomposition of Brunovsky, i.e., it decomposes the set of differential equations describing the dynamical system into two blocks of constraints: one block gathers non-linear analytical equations that do not involve differential operators, while the other block is composed of linear chains of integrators. The second part of the method, that relies on the symbolic decomposition, is numerical and based on a contractor approach. It involves a specific optimal operator for narrowing the sets of feasible solutions. This approach is shown to be efficient on a difficult problem of dynamic localization of a mobile robot, without any prior knowledge about its states.

Index Terms-Brunovsky, non-linear system, interval analysis, localization, constraint programming, differential equations.

I. INTRODUCTION

This paper introduces a new set-membership method for estimating the set of feasible trajectories of a dynamical system, without any prior knowledge about the states x(t) ∈ R n , but considering a discrete set of non-linear state observations y(t i ) ∈ R, such that:

ẋ(t) = f (x(t), u(t)) (1a) y(t i ) = g(x(t i )) (1b)
where f : R n × R m → R n is a non-linear function depicting the evolution based on an input vector u(t) ∈ R m . The observation function g : R n → R, possibly non-linear, is assumed to be scalar without loss of generality as the method is readily scalable to the vector case. The t i , i ∈ N, are observation times and the y(t i ) are the related outputs. When the states are completely unknown, conventional methods such as Kalman filters run into difficulties, as it is difficult to find a linearization point, or to perform prediction steps. Particle filters will also miss a starting point, or will employ algorithms with high complexity without ensuring a reliable convergence. In contrast, the use of set-membership approaches avoids the need for linearization and is more suited to large uncertainties by not removing consistent solutions. They are said reliable, as they will always enclose feasible solutions inside guaranteed sets. In these approaches, all errors are bounded [START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF], [START_REF] Walter | Estimation of the parameter uncertainty resulting from bounded-error data[END_REF] and an arithmetic such as interval analysis [START_REF]Applications of Interval Computations[END_REF] allows to tackle rigorously the non-linearities of the system [START_REF] Meizel | Initial localization by set inversion[END_REF], [START_REF] Desrochers | Set-membership approach to the kidnapped robot problem[END_REF].

In this context, several tools exist for estimating the set of feasible trajectories of System (1) in a reliable way. Firstly, guaranteed integration methods [START_REF] Nedialkov | VNODE-LP, A Validated Solver for Initial Value Problems in Ordinary Differential Equations[END_REF], [START_REF] Kapela | CAPD -a rigorous toolbox for Computer Assisted Proofs in Dynamics[END_REF], [START_REF] Revol | Taylor models and floating-point arithmetic: proof that arithmetic operations are validated in COSY[END_REF], [START_REF] Chapoutot | Dynibex[END_REF] can be used for dealing with the non-linear differential Eq. (1a). This amounts to solving so-called Initial Value Problems (IVPs), considering some prior state vector x 0 , possibly enclosed in a small box [x 0 ]. The outputs of these algorithms are guaranteed sets of feasible states. For many realistic applications, the tools available in the state-of-the-art are comfortable with the encountered non-linear equations. They will however badly behave in case of unknown initial condition, i.e. [x 0 ] possibly large or unbounded, and thus are not suited for the issue addressed in this paper. Secondly, the state observations y(t i ) can be processed together with their uncertainties by means of constraint propagation tools such as contractors, that are operators for safely narrowing the bounds of the state space according to a given constraint (e.g. an equation), without removing feasible solutions from the sets. However, if one considers System (1) with unbounded prior sets for enclosing the x(t) together with only few non-linear observations, the use of state-of-the-art contractor methods leads to large sets that are not relevant for state estimation purposes. In order to overcome this problem, some branch-and-prune techniques such as shaving methods [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] can be applied in order to reduce the state sets by performing bisections. The counterpart is obviously the increasing complexity of these algorithms. As a consequence, the current existing tools are not sufficient for addressing the considered problem both in a reliable and an efficient way.

Our contribution relies on the decomposition of the nonlinear dynamical system into a Brunovsky canonical form [START_REF] Brunovsky | A classification of linear controllable systems[END_REF], [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF]. The obtained decomposition is possible for flat systems which constitutes a huge class of vehicles [START_REF] Martin | Any (controllable) driftless system with m inputs and m + 2 states is flat[END_REF] such as cars [START_REF] Sekhavat | The cycab robot is a differentially flat system[END_REF], mobiles with trailers [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and applications[END_REF], sailboats [START_REF] Jaulin | Combining interval analysis with flatness theory for state estimation of sailboat robots[END_REF], etc. We will show that the obtained Brunovsky form makes it possible to solve our localization problem efficiently and in a reliable way. Besides, as this method does not involve classical interval integration tools [START_REF] Tucker | A Rigorous ODE Solver and Smale's 14th Problem[END_REF], [START_REF] Goubault | Inner approximated reachability analysis[END_REF], [START_REF] Goldsztejn | Reliable inner approximation of the solution set to initial value problems with uncertain initial value[END_REF], [START_REF] Dit Sandretto | Validated explicit and implicit Runge-Kutta methods[END_REF], [START_REF] Wilczak | Cr-lohner algorithm[END_REF], [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear hybrid systems using interval constraint propagation techniques[END_REF], the contribution of this paper brings a new approach for dealing with guaranteed integration. The paper is organized as follows. Sections II and III recall the principles of flat systems and Brunovsky decomposition. Section IV outlines set-membership methods involving contractors. It introduces the notion of Contractor Circuit which allows a better understanding of the contraction process.

Section V details a new contractor for dealing with chains of integrators, resulting from the Brunovsky decomposition. Finally, Section VI shows how to apply the proposed method for solving a state estimation problem such as the localization of a mobile robot. Convincing results are provided in this section on a reproducible robotic example. Section VII concludes the paper.

II. FLAT SYSTEMS Consider Eq. (1a) involved in the following system:

ẋ = f (x, u) z = h(x), (2)
where z ∈ R m is the output vector used with a control point of view, i.e., z is not necessarily measured. Both f and h are assumed to be smooth. The dimension of u and z are assumed to be both equal to m. The system is said to be flat with the flat output z [START_REF] Diop | Nonlinear observability, identifiability and persistent trajectories[END_REF], [START_REF] Gauthier | Deterministic observation theory and applications[END_REF], if there exists two continuous functions φ and ψ and integers κ 1 , . . . , κ m such that

   x = φ z 1 , ż1 , . . . , z (κ1-1) 1
, . . . , z m , żm , . . . , z

(κm-1) m u = ψ z 1 , ż1 , . . . , z (κ1) 1
, . . . , z m , żm , . . . , z

(3) The integers κ i , called the Kronecker indexes, are unique and should be such that [START_REF] Kalman | Kronecker invariants and feedback[END_REF]. A system satisfying Eq. ( 3) is observable: the function φ gives us the unique state vector consistent with the flat outputs and their derivatives. In practice, the functions φ and ψ involved in Eq. ( 3) can be obtained by following these steps:

κ 1 + • • • + κ m = n,
• The derivation step (see [START_REF] Isidori | Nonlinear Control Systems: An Introduction[END_REF]) computes symbolically z 1 , ż1 , . . . , z

, . . . . . . , z m , żm , . . . , z

(κm) m
as functions of x and u, using Eq. [START_REF] Walter | Estimation of the parameter uncertainty resulting from bounded-error data[END_REF]. We obtain an expression of the form

     z 1 ż1 . . . z (κm) m      = λ x u . (4)
This can be done automatically without any difficulty using Lie derivatives. It suffices to take all m equations z j = h j (x) and to compute symbolically its first, second, . . . κ j th derivatives with respect to t. At each step, the ẋi are replaced by f i (x, u).

• The inversion step has to be performed to obtain an expression of the form of Eq. (3) from the inversion of λ. This operation is difficult to obtain except for simple systems. The inversion step will not be required in our approach as it will be implicitly processed by our intervalbased resolution.

Example 1. Consider the system

   ẋ1 = x 1 + x 2 ẋ2 = x 2 2 + u z = x 1 .
(5)

For the derivation step, we compute z, ż, z, . . . with respect to x and u until u occurs. We get

   z = x 1 ż = ẋ1 = x 1 + x 2 z = ẋ1 + ẋ2 = x 1 + x 2 + x 2 2 + u. (6)
Since we had to derive twice, we conclude that the Kronecker index is κ = 2 which corresponds to the dimension of x = (x 1 , x 2 ) . As a consequence, the output z is flat. From [START_REF] Nedialkov | VNODE-LP, A Validated Solver for Initial Value Problems in Ordinary Differential Equations[END_REF] we read

λ x u =   x 1 x 1 + x 2 x 1 + x 2 + x 2 2 + u   . (7)
For the inversion step, we have to isolate x and u in order to obtain an expression with respect to z, ż, z. The inversion is easy for this example:

   x 1 = z x 2 = ż -x 1 = ż -z u = z -x 1 + x 2 + x 2 2 = z -ż -( ż -z) 2 . ( 8
)
As a consequence,

   φ (z, ż) = z ż -z ψ (z, ż, z) = z -ż -( ż -z) 2 . (9)
We can therefore consider that we have a procedure to find the functions φ and ψ as soon as we have a flat output z.

Unfortunately, to our knowledge, there is no general method to find a flat output for a given non-linear system [START_REF] Rouchon | Necessary condition and genericity of dynamic feedback linearization[END_REF]. In practice, we can use the physical intuition we have on the system [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and applications[END_REF] or tools coming from differential algebra such as the approach based on finitely generated differential ideal [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF].

III. BRUNOVSKY DECOMPOSITION

The differential flat system ẋ(t) = f (x(t), u(t)) with flat outputs z 1 , . . . , z m and sensor outputs y admits the following Brunovsky decomposition:

                                    z 1 ż1 . . . z (κm) m      = λ x u z (κ1) 1 → • • • → ż1 → z 1 . . . z (κm) m → • • • → żm → z m y = g (x) (10)
This decomposition clearly makes appear the Brunovsky normal form through the chain of integrators. It is illustrated by Fig.

(1).

Example 2. Let us consider the system described by the following equations:

       (i) ẋ1 = x 4 cos(x 3 ) (ii) ẋ2 = x 4 sin(x 3 ) (iii) ẋ3 = u 1 (iv) ẋ4 = u 2 (11) • • • • • • • • • λ z (κ m ) m z (κ 1 ) 1 z 1 z m g x y u
Fig. 1. Non-linear dynamical Systems ( 1) and ( 2) transformed into a Brunovsky canonical form. All the differential relations are gathered in chains of integrators. The rest of the system is expressed by means of non-differential functions λ and g.

This system corresponds to a wheeled robot moving on a plane. (x 1 , x 2 ) is its location, x 3 its heading, x 4 its speed and u ∈ R 2 is the control vector: acceleration and rotation speed. This system is known to be flat [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and applications[END_REF], if we take as flat output z = (x 1 , x 2 ) . To check this, we compute the Brunovsky form, applying the method of feedback linearization as explained in [START_REF] Isidori | Nonlinear Control Systems: An Introduction[END_REF]. We get that κ 1 = κ 2 = 2, and thus we have to differentiate z twice. From Eq. ( 11)-(i) and ( 11)-(ii), we obtain the following Brunovsky decomposition in which static (i.e. non-differential) constraints (I) are separated from chains of integrators (II): 12) appears as a differential equation, the expressions of each of the six variables z 1 , z 2 , . . . , z2 are considered independently and under the form of pure algebraic constraints, so that all the differential part is delegated to Block (II) in which these variables are link to each other. As expected, we have κ 1 + κ 2 = 4, which is consistent with the dimension of the system. Note that the Brunovsky Form [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF] and the state Eq.s [START_REF] Brunovsky | A classification of linear controllable systems[END_REF] are equivalent. This section has provided an illustration of how differential equations describing a dynamical system can be separated into two blocks of constraints. The first block (I) corresponds to non-linear equations with no differential operator inside, while the second block (II) is only made of chains of integrators. The following section shows how these blocks can be dealt with using contractor techniques.

(I)                                  z 1 z 2 ż1 ż2 z1 z2         =         x 1 x 2 x 4 cos(x 3 ) x 4 sin(x 3 ) u 2 cos(x 3 ) -u 1 x 4 sin(x 3 ) u 2 sin(x 3 ) + u 1 x 4 cos(x 3 )         λ(x,u) (II)    z1 → ż1 → z 1 z2 → ż2 → z 2 (12) Although Block (I) in (

IV. SET-MEMBERSHIP METHOD

The previous decomposition provides a formulation for a dynamical system which will be shown to be adapted to setmembership techniques. The relations from the decomposition are sometimes called constraints. In a set-membership context, we can formalize a state estimation problem as a Constraint Network (CN), see for instance [START_REF] Vincke | Experimental comparison of bounded-error state estimation and constraints propagation[END_REF], [START_REF] Drevelle | Reliable positioning domain computation for urban navigation[END_REF]. Interval operators, called contractors, are used to contract the set of feasible values of the state variables of the system.

A. Contractors

Contractors are algorithms based on interval analysis, that allow reliable computations despite non-linearities, large uncertainties or singularities. This arithmetic ensures to never remove valid solutions from the domains of the variables, which is critical for the safety of complex systems.

When using the framework of contractors, state estimation consists in contracting sets defined beforehand. Several contractors can be called on a set until a fixed point is reached, i.e. when the set cannot be reduced anymore. In practice, a set represents the knowledge we have about a vector or a trajectory. When nothing is known beforehand, as it is the case for the state vectors x in the considered problem, then unbounded sets are defined, such as infinite intervals or boxes [x] = [-∞, ∞] n . The same approach applies for trajectories by using tubes, i.e. intervals of trajectories denoted by [x](•), which allows to process temporal or differential equations elegantly by processing one single object. Fig. 2 illustrates a scalar tube enclosing a feasible trajectory. collection of contractors exists in the literature (see for instance [START_REF] Desrochers | A minimal contractor for the polar equation; application to robot localization[END_REF], [START_REF] Rohou | Reliable non-linear state estimation involving time uncertainties[END_REF], [START_REF] Carbonnel | Qintersection algorithms for constraint-based robust parameter estimation[END_REF]) and can be found in open-source libraries such as Codac [START_REF] Rohou | The Codac library -Constraintprogramming for robotics[END_REF]. Each of these tackles a given equation and provides an algorithm ready to be applied on sets. When a set is provided as input of such algorithm, it is possibly contracted and thus some infeasible values are removed. For instance, the equation f (x) = 0 can be addressed by a contractor C HC4revise provided in [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] and applied on some box

t [x] t0 tf [x]([t1]) [t1] [x](t2) t2 x * ( •) x -(•) x + (•)
[x]. The obtained contracted box C HC4revise ([x]
) is a subset of [x] but still contains all feasible vectors consistent with f (x) = 0. This consistency property of contractors allows to apply them as much as is necessary. Several contractors can be combined in order to tackle more complicated problems. From Tarski's theorem [START_REF] Tarski | A lattice-theoretical fixpoint theorem and its applications[END_REF] and since all contractors we use are inclusion monotonic, the resulting propagation method always converges.

B. Contractor Circuits

Successive contractions allow to converge towards accurate outer approximations of sets of feasible solutions for the considered problem. This is achievable by combining several contractors on common sets. Usually, a decomposition of the problem is performed in order to bring out equations associated with already available contractors. For instance, the Brunovsky Form (10) already provides a decomposition for which contractors are at hand. Graphically, the contractor resolution amounts to the Contractor Circuit depicted in Fig. 3. 3. So-called Brunovsky Circuit associated with the Brunovsky decomposition of the robotic Eq. ( 12).

[x](•) [u](•) [z 1 ](•) [ ż1 ](•) [z 1 ](•) [z 2 ](•) [ ż2 ](•) [z 2 ](•) C C λ C Fig.
The concept of Contractor Circuit, introduced in this paper, allows a better understanding of the use of contractors for state estimation. A Contractor Circuit is a Boolean matrix that links contractors to sets. In this representation, entries of the matrix are represented by the small black dots at the intersection of horizontal and vertical lines. At the bottom of the figure, each rectangular block corresponds to a contractor involved in the circuit. Horizontal lines are associated with sets, i.e. tubes in this example. The red horizontal lines correspond to the inputs and the blue ones to the output variables. The black vertical lines represent the possible inputs of the contractor algorithms.

The initial bounds of the tubes correspond to some prior knowledge we have on the values, or to some measurements that have been performed. The latter can be seen as a restriction on the sets: the measurement is directly encoded in the set. Tubes are therefore initialized to [-∞, ∞], ∀t if there is not prior information available for the related variables, as it is the case for the state trajectory x(•). The contractors involved in the circuit are then called until no more contraction can be observed. Note that the order in which the contractors are called will not affect the final results [START_REF] Apt | Principles of constraint programming[END_REF], because all the operators are called indefinitely until a fixed point is reached. The use of a Boolean matrix can be avoided but the call of the contractors and the connections between tubes should then be done by hand inside the program. Using the connection matrix allows us to enclose a part of the usual code inside the data structure of the problem. As a consequence, it helps to define a high level language for computing with contractors.

The Circuit depicted in Fig. 3 involves the contractor C λ related to the non-linear equation (z, ż, z) = λ(x, u). The literature already provides algorithms for addressing this relation, see for instance the HC4revise contractor introduced in [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF]. In addition, the circuit involves two C operators that are the contractors for dealing with integrator chains. Because the dimension of z is equal to 2, two contractor chains are involved in this circuit. C cannot be directly found in the literature, but is required for Brunovsky decompositions. The following section proposes an algorithm for C .

V. THE INTEGRATOR CHAIN CONTRACTOR

In order to process the Brunovsky decomposition in a setmembership context, a dedicated integrator chain contractor, denoted by C , has to be provided. As explained in Section III, the integrator chain encompasses all differential links of System (1) in one single relation given by:

z (κ) → • • • → ż → z ( 13
)
where z is a scalar output as defined in System (2). For vector outputs z, Relation (13) is involved several times according to the dimension of z. This integrator chain allows to accurately propagate information from one signal through its primitives and derivatives. Note that this relation could be broken down into several elementary integrators as in [START_REF] Sekhavat | The cycab robot is a differentially flat system[END_REF], for which contractors already exist, see for instance [START_REF] Rohou | Guaranteed computation of robot trajectories[END_REF].

z (κ) → z (κ-1) , . . . , z → ż, . . . , ż → z. (14)
However, in a set-membership context, this leads to wellknown wrapping effects due to the enclosure of each intermediate result into non-optimal sets. The use of one single operator allows optimal computations.

A. Linear state estimator

The chain (13) can be expressed under the form of a specific linear state estimation for which the state is defined by a set z of primitives/derivatives of one only signal 1 and for which the observations are provided by restrictions on these primitives/derivatives. More precisely, the integrator chain constraint involving the signals z (0) , z (1) , . . . , z (κ) , w and defined as:

w → z (κ) → • • • → ż → z ( 15
)
can be cast into the following linear system corresponding to the Brunovsky form:

ż(t) =        0 1 0 • • • 0 0 1 . . . . . . . . . 0 0 0 • • • 1 0 0 0 • • • 0        A z(t) +        0 0 . . . 1        B w(t), (16)
where w(•) is known to be inside a tube [w](•). An integrator chain contractor 1) ](•), . . . , [z (κ) ](•) from Eq. ( 16) and [w](•).

C [z](•), [w](•) has to be built to contract a tube [z](•) = [z (0) ](•), [z (
In the set-theoretical framework, due to the linearity of System ( 16), an optimal contraction is achievable by computing a set of abstract sets Z(t) over time, using the following equations:

-→ Z (t) ⊂ Ž(t) ∩ e Aδ • -→ Z (t -δ) + δe A•[0,δ] • B • [w]([t -δ, t]) ← - Z (t) ⊂ Ž(t) ∩ e -Aδ • ← - Z (t + δ)-δe -A•[0,δ] • B • [w]([t, t + δ]) Z(t) = -→ Z (t) ∩ ← - Z (t)
where • Ž(t) is the prior set of states for t;

• - → Z (t),
← -Z (t) are the sets of all states at t consistent with the past and the future, respectively; • Z(t) is the smallest set for z(t) assuming that w(t) ∈

[w](t) and z(t) ∈ Ž(t). It is the result of the contraction.

It has been shown that this observer is exact in the discrete time case, if we represent the sets Ž(t), -→ Z (t), ← -Z (t), Z(t) by polygons [START_REF] Walter | Exact recursive polyhedral description of the feasible parameter set for bounded-error models[END_REF], [START_REF] Wan | Guaranteed state estimation for nonlinear discrete-time systems via indirectly implemented polytopic set computation[END_REF]. An implementation in the continuous-time case has been recently proposed in [START_REF] Rohou | Exact bounded-error continuous-time linear state estimator[END_REF]. The principle is to use a tube for representing prior sets Ž(•) and to apply the exact discrete time observer of [START_REF] Walter | Exact recursive polyhedral description of the feasible parameter set for bounded-error models[END_REF] for each slice of the tube. The resulting contractor is still minimal if we assume that we can neglect the discretization errors, otherwise, an unavoidable small pessimism exists, due to the implementation of tubes.

Therefore, state observations are processed as restrictions on the Ž(t) prior sets. In our approach, these restrictions directly come from the tubes, under the form of boxes, since our tubes are implemented as sets of boxed slices. The contraction results into the intersection between the integration of the prior sets Ž(tδ) into -→ Z (t) (from Eq. ( 16)) and the prior sets Ž(t). It always corresponds to the intersection of a polygon and a box, which can be computed accurately. Fig. 4 provides a simple illustration of this state estimation in the continuous case, using polygons and tubes. Considering for instance the chain w → z 2 → z 1 , prior 2d sets Ž(•) are implemented as tubes [z 1 ]×[z 2 ](•) (upper part of the figure). The lower part of the figure depicts in gray one integration step (polygon Ž(t 1 ) into polygon -→ Z (t 1 + δ)) and the intersection with the prior set at time t 1 + δ: the box Ž(t 1 + δ). The resulting contracted set, represented in hatched blue, is also a polygon.

The guaranteed linear integration transforms each vector into a box and each box or polygon into a polygon. The use of polygons for representing sets allows to obtain optimal results

z 1 (t 1 ) z 2 (t 1 ) z 1 (t 1 + δ) z 2 (t 1 + δ)
Eq. ( 16) with minimal enclosures. However, it must be highlighted that this optimality comes with a strong complexity of the algorithms, since the number of vertices increases dramatically with time. Non-minimal enclosures could be used instead, such as intervals [START_REF] Xiong | Fault detection using interval kalman filtering enhanced by constraint propagation[END_REF], [START_REF] Mustafa | Guaranteed SLAM; An Interval Approach[END_REF], zonotopes [START_REF] Combastel | A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes[END_REF] or parallelotopes [START_REF] Chisci | Recursive state bounding by parallelotopes[END_REF], allowing faster computations. In practice, these enclosures would be more suited for state estimation, since the optimality of this differential contractor is not always relevant with respect to the number of observations of the system. A first study has recently led to the publication of [START_REF] Rauh | An ellipsoidal predictor-corrector state estimation scheme for linear continuous-time systems with bounded parameters and bounded measurement errors[END_REF] in which an ellipsoidal contractor is provided for ż = Az + Bw with higher performances (faster computations, implementation for any dimension) than the polygonal approach of [START_REF] Rohou | Exact bounded-error continuous-time linear state estimator[END_REF], while maintaining a quite accurate envelope of the states. Its application in the context of this paper remains to be studied, but is encouraging for real time applications.

t z 2 (•) t z 1 (•) δ δ t1 t1 + δ t1 t1 + δ ob se rv er a [z1](t1) × [z2](t1) Ž1(•) Ž2(•) Ž(t1) - → Z (t1 + δ) Z(t1 + δ) Ž(t1 + δ)
As a conclusion, a contractor dealing with linear systems expressed as ż = Az + Bw can be configured in order to deal with a chain of integrators. This allows to contract a list of tubes directly linked by their derivative and primitive relations. In the remaining of this paper, the contractor provided in the publication [START_REF] Rohou | Exact bounded-error continuous-time linear state estimator[END_REF] is used for the C involved in the Brunovsky circuit.

VI. APPLICATION TO LOCALIZATION

The set-membership Brunovsky decomposition proposed in this paper is now illustrated on a difficult academic robotic problem.

A. Problem statement

We consider the localization problem of a robot moving in an environment where few landmarks are available [START_REF] Leonard | Mobile robot localization by tracking geometric beacons[END_REF], and formalized by System (1). This state estimation is difficult if all state vectors x(t) ∈ R 4 , describing the positions (x 1 , x 2 ) , headings x 3 and speeds x 4 of the vehicle are unknown, i.e. when initial conditions or direct measurements of the states are not at hand. Only the control u(t) ∈ R2 of the system (linear acceleration and rotation rate) is known, with some bounded uncertainties, together with some interval range-only observations [y](t i ) made between the moving vehicle and two landmarks, for which the positions m a and m b are known.

The goal of this problem is to perform a full state estimation by computing an approximation of the four components of x(t) along time, including the headings x 3 (t). Such state estimation is difficult due to the absence of knowledge on the initial condition x(0) as well as any other state, and because very few asynchronous non-linear observations are available. In particular, it has been shown in [START_REF] Wang | Dynamic icsp graph optimization approach for car-like robot localization in outdoor environments[END_REF] and [START_REF] Colle | Mobile robot localization by multiangulation using set inversion[END_REF], that the problem is difficult to solve if the heading of the robot is not measured. As it has been said before, classical methods such as Bayesian approaches or particle filters badly behave in presence of few observations and are not suited to deal with large uncertainties because of inaccurate linearizations or convergence issues. On the other hand, set-membership approaches are more dedicated to limited observations and non-linear contexts, but very few studies have been done on this application [START_REF] Vincke | Experimental comparison of bounded-error state estimation and constraints propagation[END_REF]. One should note the work of [START_REF] Wang | Guaranteed localization and mapping for autonomous vehicles[END_REF] based on interval methods with a pruning component [START_REF] Janssen | Optimal Pruning in Parametric Differential Equations[END_REF], that we want to avoid. Indeed in [START_REF] Wang | Guaranteed localization and mapping for autonomous vehicles[END_REF], the approximation of heading values is obtained by bisecting the heading space. Therefore, a large number of bisections is needed in order to obtain a significant approximation [START_REF] Kieffer | Nonlinear identification based on unreliable priors and data, with application to robot localization[END_REF], which ruins the performances of interval methods [START_REF] Guyonneau | A visibility information for multi-robot localization[END_REF], [START_REF] Porta | A Space Decomposition Method for Path Planning of Loop Linkages[END_REF].

B. Reproducible test case

Let us consider the following reproducible example. The proposed simulation is based on analytical expressions and simple data in order to encourage future comparisons with the method provided in this paper. A wheeled robot described by Eq. ( 11) follows a planar trajectory given by a Lissajous equation. The unknown states are given by:

x(t) =    
10 cos(t) 5 sin(2t) atan2 (10 cos(2t), -10 sin(t)) (-10 sin(t)) 2 + (10 cos(2t))

2     . (17)
Eq. ( 17) cannot be used during the resolution. It is only given to validate the state estimation, and also to generate data (exteroceptive measurements and inputs). From Eq. ( 17), we obtain the following inputs:

u(t) =   2 sin(t) sin(2t)+cos(t) cos(2t) sin 2 (t)+cos 2 (2t) 10 cos(t)•sin(t)-20 cos(2t)•sin(2t) √ sin 2 (t)+cos 2 (2t)   . ( 18
)
The inputs are supposed accurately known. In our method, they are actually numerically enclosed in a thin tube [u](•) that is not degenerated because of the sliced implementation of the tubes of the used library [START_REF] Rohou | The Codac library -Constraintprogramming for robotics[END_REF]. Finally, we assume that the robot performs distance measurements between its position and one of the landmarks m a or m b , respectively located at (-5, 6) and (0, -4). The observation equation (1b) is therefore the distance function:

y j (t i ) = x 1 (t i ) -m j 1 2 + x 2 (t i ) -m j 2 2 . ( 19
)
The beacons emit asynchronously range-only signals received by the robot. These bounded measurements are given in Table I. The simulation is run for t ∈ [0, 3].

C. Interval Brunovsky localization: methodology

Our interval state estimation approach consists in the following steps:

• First, a symbolic decomposition of Eq. ( 11) is made in order to obtain a Brunovsky form. The decomposition makes appear two blocks of differential and nondifferential relations, as in [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF].

• Second, sets of feasible values are defined in order to enclose the involved variables x(•), u(•), and the [y](t i )'s.

For the states, we define an unbounded tube

[x](•) such that ∀t, [x](t) = [-∞, ∞] 4 . The tube [u](•)
is thin and centered on Eq. ( 18). This means that the analytical knowledge of u(•) is not required by our method 2 . Finally, interval measurements [y](t i ) are given by Table I.

• The third step consists in defining contractors for the equations resulting from the first step. This amounts to build the Brunovsky Circuit of Fig. 3. The non-differential and non-linear function λ(x, u) is given by ( 12) and a contractor C λ can be easily built from this analytical expression [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF]. Then, the two C of Fig. 3 are contractors for dealing with linear observers, such as those provided in [START_REF] Rohou | Exact bounded-error continuous-time linear state estimator[END_REF] or [START_REF] Rauh | An ellipsoidal predictor-corrector state estimation scheme for linear continuous-time systems with bounded parameters and bounded measurement errors[END_REF], and configured according to System [START_REF] Jaulin | Combining interval analysis with flatness theory for state estimation of sailboat robots[END_REF]. We add another contractor C g for dealing with Eq. ( 19), also built from an analytical expression. This last contractor links some slices [x 1,2 ](t i ) of the tube [x](•) with the measurements [y](t i ). More details are given in [START_REF] Rohou | Reliable non-linear state estimation involving time uncertainties[END_REF].

• Finally, the contractors are applied on the intervals and tubes defined in the second step. Some iterative resolution is performed by calling all the contractors until no more contraction is obtained. It should be recalled that the order of calls does not change the quality of the results, only the computation time may differ. Once this fixed point has been reached, the estimated (feasible) states are enclosed in the contracted tube [x](•).

D. Results

The obtained results for the position tube [x 1,2 ](•) is illustrated in blue in Fig. [START_REF] Desrochers | Set-membership approach to the kidnapped robot problem[END_REF]. It encloses the unknown true trajectory given by Eq. ( 17) that is white painted. As one can see, despite the absence of prior knowledge on the states, the obtained envelope is thin thanks to the use of efficient contractors. The two landmarks are pictured by orange boxes and the four range-only measurements are drawn by gray lines. The simulation, run in 36 seconds, also provides a contraction of the velocity tube [x 4 ](•) and the heading tube [x 3 ](•), as illustrated by Fig. 6. The example has been processed by operators delivered in the Codac library [START_REF] Rohou | The Codac library -Constraintprogramming for robotics[END_REF], in which this example is freely available. The Reader may find more material related to this work on http://codac.io/brunovsky.

VII. CONCLUSIONS

This paper has presented a new approach for state estimation of non-linear dynamical systems, that combines some non-linear state estimation techniques [START_REF] Fliess | Control via state estimations of some nonlinear systems[END_REF] based on flatness [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and applications[END_REF] with interval contractors. Flatness allows us to build a Brunovsky decomposition made with a set of static nondifferential and non-linear relations, together with chains of integrators. A contractor based approach is then able to take advantage of this decomposition in order to enclose, in a guaranteed way, the unknown states into an accurate tube.

In our set-membership approach, all the differential relations are processed by one single operator, able to deal with a chain of integrators, and providing optimal results when using polygons as wrappers. The resulting state estimator has several advantages over classical approaches:

• The state estimator is reliable with respect to nonlinearities. Thanks to interval analysis, it is able to deal with non-linear (or non-differentiable and even noncontinuous) state equations, without linearizing (as done by the extended Kalman filter [START_REF] Thrun | Probabilistic Robotics[END_REF]) or approximating them.

• The state estimator does not require the interval integration of differential equations. Such integrations are needed by all other interval state estimation methods [START_REF] Wilczak | Cr-lohner algorithm[END_REF], [START_REF] Revol | Taylor models and floating-point arithmetic: proof that arithmetic operations are validated in COSY[END_REF], [START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear differential equations[END_REF], [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], which makes them inefficient for large uncertainties.

• The state estimator can be used for real-time applications. Our approach has been applied to the range-only localization of a robot, without any measurements of its heading or speed, nor any knowledge about its positions. This state estimation problem, that is known to be difficult to address, has been accurately and reliably solved using the Brunovsky decomposition and the specific operators described in this paper, without any bisection technique or prior knowledge on the states.

Fig. 2 .

 2 Fig. 2. A tube [x](•) defined on the temporal domain [t 0 , t f ] as an interval of two signals [x -(•), x + (•)] and enclosing an unknown trajectory x * (•). The thinner the tube, the better the approximation of x * (•).

Fig. 4 .

 4 Fig. 4. Example of integrator chain contractor over one time step.

Fig. 5 .

 5 Fig. 5. Set [x](•) of feasible states projected in two dimensions. The unknown planar trajectory remains enclosed in the tube, which illustrates the guarantee of the computations.

Fig. 6 .

 6 Fig. 6. The tube [x 3 ](•) of feasible headings. The actual but unknown truth, given by x * 3 (t) = atan2 10 cos(2t), -10 sin(t) , is plotted in white and guaranteed to be enclosed in the computed [x 3 ](•).

TABLE I SET

 I OF FOUR BOUNDED MEASUREMENTS (t i , [y](t i )).

	t i	[y a ](t i )	t i	[y b ](t i )
	0.75	[12.333,12.383]	1.50	[4.733,4.783]
	2.25	[10.938,10.988]	3.00	[10.211,10.261]

The notation z is chosen voluntarily as it will refer, in the application of the contractor related to this state estimator, to the output of System (2).

Because our tubes are implemented as lists of boxes, some pessimism related to wrapping effects is introduced. The time-discretization of our tubes is set to δ = 0.003.

ACKNOWLEGMENT

This work was supported by the French Agence Nationale de la Recherche (ANR) [grant number ANR-16-CE33-0024].