C V Mai

B Sudret

Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping

Keywords:

Polynomial chaos expansions (PCE) have proven efficiency in a number of fields for propagating parametric uncertainties through computational models of complex systems, namely structural and fluid mechanics, chemical reactions and electromagnetism, etc. For problems involving oscillatory, time-dependent output quantities of interest, it is well-known that reasonable accuracy of PCE-based approaches is difficult to reach in the long term. In this paper, we propose a fully non-intrusive approach based on stochastic time warping to address this issue: each realization (trajectory) of the model response is first rescaled to its own time scale so as to put all sampled trajectories in phase in a common virtual time line. Principal component analysis is introduced to compress the information contained in these transformed trajectories and sparse PCE representations using least angle regression are finally used to approximate the components.

The approach shows remarkably small prediction error for particular trajectories as well as for second-order statistics of the latter. It is illustrated on different benchmark problems well known in the literature on time-dependent PCE problems, ranging from rigid body dynamics, chemical reactions to forced oscillations of a non linear system.

Introduction

In modern engineering, it is of utmost importance to investigate the significant effects of uncertainties when considering the behaviour of complex systems. These uncertainties may arise from 1 environmental factors (e.g. excitations, boundary conditions) or inherent sources (e.g. natural variability of the materials) and are usually represented by random variables. In this context, the framework of uncertainty quantification was introduced, of which a major component is the propagation of uncertainty from the input parameters through the system to the output quantities of interest. The outcomes of uncertainty propagation (e.g. statistical, reliability and sensitivity measures) allow a better understanding of the system and are critical in decision making.

So far Monte Carlo simulation (MCS) is universally used for solving uncertainty propagation problems. The idea behind MCS is to perform the simulation a sufficiently large number of times by varying input parameters such that the average of the response quantity of interest converges to the expected value according to the law of large numbers. However, the use of MCS is hindered by the fact that a large number of simulations is not affordable in many practical problems (e.g. when each evaluation of the computational model is time-and memory-consuming).

To overcome this issue, spectral methods have been used in the last two decades as an alternative approach to traditional MCS. The spectral approach consists in representing the response quantity of interest in a space spanned by well-defined basis functions. Among a wide variety of basis functions that have been investigated, polynomial functions have shown particular effectiveness [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]Le Maître and Knio, 2010;[START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]. The spectral approach that uses polynomial chaos functions as a basis is simply named polynomial chaos expansions (PCEs) in the sequel.

In practice, PCEs are widely used as an approximate model to substitute a computationally expensive model for uncertainty propagation. They can be used in an either intrusive or nonintrusive setup. The former requires knowledge of the mathematical equations describing the considered system. One has to interfere with the original set of equations, reformulate it and then solve the reformulated system to compute the PCE coefficients. In contrast, the latter does not necessitate any prior knowledge of the governing equations. It considers the deterministic computational model as a blackbox and only requires to define an experimental design, i.e. a set of input and corresponding output values. In several studies, PCEs have shown great efficiency compared to the traditional uncertainty propagation approach with MCS, see e.g. Dossantos-Uzarralde and Guittet (2008); [START_REF] Rajabi | Polynomial chaos expansions for uncertainty propagation and moment independent sensitivity analysis of seawater intrusion simulations[END_REF].

PCEs, however, face challenges when used for dynamical systems that are encountered in the fields of structural and fluid dynamics or in chemical engineering (Wan and Karniadakis, 2006a;Beran et al., 2006a;[START_REF] Ghosh | Applicability of the spectral stochastic finite element method in time-dependent uncertain problems[END_REF]Le Maître et al., 2010). In these cases, the governing equations are a system of ordinary differential equations with random parameters.

First, the response as a function of time is no longer a scalar quantity but may be cast as a vector after proper time discretization. Applying PCEs at each time instant might require large computational resources. To reduce the computational cost, [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF] used principal component analysis to capture the main stochastic features of the vector-valued response quantities by means of a small number of variables which can be represented by PCEs.

The greatest challenge is the decrease in time of the accuracy of the PCE model as reported in numerous publications (Wan and Karniadakis, 2006a;Beran et al., 2006a;[START_REF] Ghosh | Applicability of the spectral stochastic finite element method in time-dependent uncertain problems[END_REF]Le Maître et al., 2010;Le Maître and Knio, 2010;[START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] though. The features of the accuracy degeneration, i.e. its onset (the instant at which PCEs start being insufficiently accurate) or its rate (how fast the accuracy decreases), depend on the considered problem.

The cause of the decaying accuracy of PCEs in dynamics can be classified into an approachrelated cause and an inherent cause. The approach-related cause refers uniquely to intrusive techniques. In fact, the latter solves a system of reformulated ordinary differential equations which are derived from the original system of equations by substituting PCE for the quantity of interest. At any given instant, the PCE is truncated after P terms, thus introducing a truncation error. The latter is accumulated in time, therefore the results deteriorates [START_REF] Ghosh | Applicability of the spectral stochastic finite element method in time-dependent uncertain problems[END_REF]. By means of the non-intrusive approach, one can avoid this source of error since the responses at different instants can be examined "independently", which prevents the accumulation of error at later instants provided the deterministic solver is equally accurate whatever the realization of the input parameters. The inherent cause refers to the fact that the problem itself demonstrates increasing complexity as time evolves, as shown through examples in [START_REF] Pettit | Spectral and multiresolution Wiener expansions of oscillatory stochastic processes[END_REF]; [START_REF] Ghosh | Applicability of the spectral stochastic finite element method in time-dependent uncertain problems[END_REF]. The growing complexity makes it increasingly hard for PCE to capture the behaviour of the system.

The growth in time of the inherent complexity of the problem is characterized by the increasingly complicated relationship between the output quantity and the input parameters, exhibiting important non-linearity, abrupt changes and possibly discontinuities (see e.g. [START_REF] Witteveen | Simplex stochastic collocation with eno-type stencil selection for robust uncertainty quantification[END_REF]; [START_REF] Desai | Uncertainty quantification of a nonlinear aeroelastic system using polynomial chaos expansion with constant phase interpolation[END_REF]). It may be related to the difference in terms of frequency and phase content of the various response time-series obtained with distinct values of the uncertain input parameters (Le Maître et al., 2010;[START_REF] Witteveen | An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase[END_REF]. These discrepancies tend to be more and more severe when time passes. In other words, trajectories tend to be similar at early instants and less and less in phase in the long term (Wan andKarniadakis, 2005, 2006a).

Existing approaches available in the literature for dealing with dynamical systems subject to uncertainties are summarized hereafter. For the sake of clarity, the approaches are classified as shown in Figure 1.

As explained previously, the model complexity increases in time, thus the PCE degree is required to increase accordingly in order to achieve an acceptable level of accuracy. In high dimensional problems, the use of high-order PCEs leads to the exponential increase of the size of the basis. However, the amount of available information, i.e. the size of the experimental design used for fitting PCEs, usually remains limited. This approach will not be effective after a certain period of time because of this so-called curse of dimensionality.

To alleviate this issue, [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] introduced adaptive sparse PCEs that allows one to take advantage of the sparsity in the structure of the model (if this sparse structure exists), thus extending the time range where the computation of PCEs is tractable and the result is sufficiently accurate. In other words, adaptive sparse PCEs may delay the onset of the accuracy degeneration. [START_REF] Lucor | Adaptive generalized polynomial chaos for nonlinear random oscillators[END_REF] used adaptive generalized PCEs, which consists in detecting the first-order terms with the most important effects on the fluctuation of the response and then building the higher-order terms that only include the selected first-order terms. From the same perspective, Mai and Sudret (2015a) developed the hierarchical PCEs which aims at updating the set of candidate polynomials adaptively by adding selected interaction terms while selecting only the regressors with the most importance. In most papers, the proposed high-order PCE approaches consist in using assumptions to reduce the size of the high-order PCE basis or using advanced computational techniques for computing them. Wan andKarniadakis (2005, 2006a,b) proposed multi-element PCEs, in which the random space is divided into multiple subspaces in such a way that the complexity of the model in each subspace is reduced, thus requiring only low-order PCEs. [START_REF] Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF] also used multielement PCEs with a discontinuity detector in order to minimize the number of subspaces. [START_REF] Nouy | Identification of multi-modal random variables through mixtures of polynomial chaos expansions[END_REF] and [START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] approximated a multimodal random variable (i.e. the output quantity of interest) by a mixture of unimodal random variables, each modelled by PCEs. This approach might help to improve the effectiveness of PCEs in the context of dynamical systems, when the responses at late instants usually exhibit multi-modal distributions as will be shown in the current paper through numerical applications. In the above approaches, the input space is divided into subspaces according to the detected discontinuities or dissimilarities. One then builds a local PCE in each subspace and combines those PCE models to obtain a global metamodel.

Therefore these approaches can be classified as local PCEs. The use of polynomial functions in local domains, however, requires an accurate decomposition of the input space and will not be straightforward in high-dimensional problems.

From a different perspective, [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] proposed to compute time-dependent

PCEs by updating the polynomial chaos basis on-the-fly. If the approximation error is excessive at a considered time instant, the authors add to the existing set of random variables a new variable, which is the response quantity at the previous instant. This is based on the idea that a fixed set of random variables at the beginning of the process is not sufficient to model the system in the long term and thus, the set of random variables the PCEs depend on needs to be updated. This approach can be viewed as a nested PCE model, i.e. a PCE model of another PCE model. [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF] used flow map composition, which is in principle similar to time-dependent PCEs. The time-history response is composed of short-term flow maps, each modelled by PCE. The idea of constructing the basis on-the-fly was also applied by [START_REF] Cheng | A dynamically bi-orthogonal method for timedependent stochastic partial differential equations I: Derivation and algorithms[END_REF] and [START_REF] Choi | On the equivalence of dynamically orthogonal and bi-orthogonal methods: Theory and numerical simulations[END_REF], who derived intrusively a system of equations governing the evolution of the time-dependent spatial and stochastic basis. In the context of structural dynamics, Spiridonakos and Chatzi (2015a,b) proposed the combination of PCEs and autoregressive models which consists in representing the response as a function of its past values.

This approach is currently investigated with the use of sparse adaptive PCEs by [START_REF] Mai | Surrogate modelling for stochastic dynamical systems by combining NARX models and polynomial chaos expansions[END_REF].

Recently, [START_REF] Ozen | Dynamical polynomial chaos expansions and long time evolution of differential equations with random forcing[END_REF] introduced the dynamical PCEs, which is also based on the idea that the future evolution of the response depends on the present solution.

As explained earlier, the accuracy of PCEs may degenerate in time due to the time-increasing dissimilarity between the response trajectories when considering distinct values of the uncertain input parameters. To alleviate the accuracy decay, one may naturally think of increasing the similarity between the response trajectories. For this purpose, an attractive approach is to preprocess the response trajectories in order to increase the similarity between them. To this end, [START_REF] Witteveen | An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase[END_REF]; [START_REF] Desai | Uncertainty quantification of a nonlinear aeroelastic system using polynomial chaos expansion with constant phase interpolation[END_REF] As a summary, PCEs fail to represent long-term time-dependent system responses because of their inherent increasing complexity. To the authors' knowledge there is no versatile tool that helps overcome the problem in a non-intrusive setup. This paper aims at filling this gap by introducing a fully non-intrusive approach that allows efficient use of PCEs for time-dependent problems showing oscillatory behaviours. The proposed approach relies on a stochastic time warping and the subsequent rescaling of the response trajectories.

The paper is organised as follows: in Section 2, the fundamentals of PCEs for time-independent problems are recalled. We introduce so-called time-frozen PCEs that will be used for comparison. In Section 3, we propose an original non-intrusive PCE approach for uncertain dynamical systems based on stochastic time-warping. Five applications are finally considered to show the efficiency of the proposed approach.

2 Polynomial chaos expansions

Spectral representation

Let us consider the model Y = M(X) where X = (X 1 , . . . , X M ) is a M -dimensional input vector of random variables with given joint probability density function f X defined over an underlying probability space (Ω, F, P) and M : x ∈ D X ⊂ R M → R is the computational model of interest, where D X is the support of the distribution of X. Herein, we assume that the input random variables are independent, i.e. the joint probability density function (PDF) is the product of the marginal PDFs:

f X (x) = f X1 (x 1 ) . . . f XM (x M ). (1)
Assuming that the scalar output Y is a second order random variable, i.e. E Y 2 < +∞, is equivalent to require that the computational model M belongs to the Hilbert space H of square-integrable functions with respect to the inner product:

< u , v > H = DX u(x)v(x)f X (x)dx. ( 2
)
Denote by H i the Hilbert space of square-integrable functions with respect to the marginal probability measure P Xi (dx i ) = f Xi (x i )dx i . Let us equip H i with an inner product:

< u , v > Hi = DX i u(x i )v(x i )f Xi (x i )dx i , (3)
where D Xi is the support of the distribution of X i and denote by {φ i k , k ∈ N} an orthonormal basis of H i which satisfies:

< φ i k , φ i l > Hi = δ kl , (4)
in which δ kl is the Kronecker symbol, which is equal to 1 if k = l and equal to 0 otherwise.

As shown by [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF], the Hilbert space H is isomorphic to the tensor product ⊗ M i=1 H i . Thus a basis of H may be obtained by the tensor product of the univariate bases φ i k , k ∈ N , i = 1, . . . , M . As a consequence, the random variable Y = M(X) that results of the propagation of the uncertainties modelled by X through the computational model M may be cast as:

Y = α1∈N . . . αM ∈N y α1...αM φ 1 α1 (X 1 ) . . . φ M αM (X M ). ( 5
)
For the sake of simplicity, introducing multi-indices α = {α 1 , . . . , α M }, Y may be rewritten as:

Y = α∈N M y α φ α (X). (6)
where φ α (X) = M i=1 φ i αi (X i ) are the multivariate basis functions and y α are the associated deterministic coefficients.

Polynomial chaos expansions

Univariate basis functions φ i k , k ∈ N, i = 1, . . . , M may be constructed using orthonormal polynomials [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] leading to the so-called generalized polynomial chaos expansion [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]. For instance when X i is a uniform (resp. standard normal) random variable, the corresponding polynomial basis comprises orthonormal Legendre (resp. Hermite) polynomials. Then Eq. ( 6) becomes:

Y = α∈N M y α ψ α (X), (7)
in which α = (α 1 , . . . , α M ) are the multi-indices with α i , i = 1, . . . , M denoting the degree of the univariate polynomial in X i and ψ α (X) = M i=1 ψ i αi (X i ) are multivariate orthonormal polynomials obtained by the tensor product of univariate polynomials.

In practice, the use of infinite-dimensional PCEs is not tractable. One always truncates the expansion to obtain an approximate representation:

Y = α∈A y α ψ α (X) + , (8)
in which A is a truncation set and is the truncation-induced error. A classical truncation scheme consists in selecting all polynomials of total degree less than or equal to p, when the truncation set reads:

A M,p = {α ∈ N M : α 1 def = α 1 + . . . + α M p}. (9)

Computation of PC coefficients and error estimation

The computation of the coefficients {y α , α ∈ A} in Eq. ( 8) can be conducted using intrusive (i.e. Galerkin scheme) or non-intrusive approaches (e.g. projection, regression and quadrature methods). In the following, we will compute the coefficients of the expansions using the adaptive sparse PCE technique proposed by [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] which is a non-intrusive leastsquare minimization technique based on the least angle regression algorithm [START_REF] Efron | Least angle regression[END_REF].

The reader is referred to [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] for more details on this approach.

The accuracy of the representation is estimated by means of the leave-one-out (LOO) crossvalidation, which allows a fair error estimation at an affordable computational cost [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]. The principle of cross validation is to use different sets of points to build PCEs, then compute the errors with the actual model. Assume that one is given a sample set X = x (i) , i = 1, . . . , n . The computational model M is run for each point in X , resulting in the vector of output quantity values Y = y (i) , i = 1, . . . , n . Setting one point x (i) apart from X , one can build a PCE model M PC\i (•) from the remaining points X \x (i) = x (1) , . . . , x (i-1) , x (i) , . . . , x (n) . The predicted residual error at point x (i) reads:

∆ (i) def = M(x (i) ) -M PC\i (x (i) ). ( 10
)
The LOO error is defined as follows:

Err LOO = 1 n n i=1 ∆ 2 i . (11)
At first glance, one could think that evaluating the LOO error is computationally demanding since it requires n different predicted residuals, each of them obtained from a different PCE.

However, by means of algebraic derivations, one can compute Err LOO from a single PCE M PC (•) built with the full experimental design as follows [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]:

Err LOO = 1 n n i=1 M(x (i) ) -M PC (x (i) ) 1 -h i 2 , ( 12
)
where h i is the i th diagonal term of the projection matrix A A T A -1

A T and the information matrix A is defined by A ij = ψ j (x (i) ), i = 1, . . . , n, j = 1, . . . , card A , i.e. the i th row of A is the evaluation of the polynomial basis functions at the point x (i) in the ED. Note that in practice, a normalized version of the LOO error is used:

ˆ LOO = Err LOO Var [Y] , (13)
where Var [Y] is the empirical variance of the sample of outputs.

Time-frozen polynomial chaos expansions

In the context of time-dependent problems, i.e. Y (t) = M(X, t), the polynomial chaos representation of the response quantity reads:

Y (t) = α∈A y α (t)ψ α (X) + (t) (14)
in which the notation y α (t) indicates the time-dependent coefficients of PCEs. The representation of a time-dependent quantity by means of PCEs as in Eq. ( 14) is widely used in the literature, see e.g. [START_REF] Pettit | Spectral and multiresolution Wiener expansions of oscillatory stochastic processes[END_REF]; Le Maître et al. ( 2010); [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF].

At a given time instant t, the coefficients {y α (t), α ∈ A} and the accuracy of the PCEs are estimated by means of the above mentioned techniques (see Section 2.3). The metamodel of the response is computed independently at each time instant, hence the name time-frozen PCEs.

We now introduce the use of time-frozen PCEs for computing the time-dependent statistics of the response. The multivariate polynomial chaos functions are orthonormal, i.e. :

E ψ α (X) ψ β (X) def = DX ψ α (x) ψ β (x) f X (x) dx = δ αβ ∀α, β ∈ N M , (15)
in which δ αβ is the Kronecker symbol that is equal to 1 if α = β and equal to 0 otherwise. Indeed,

each multivariate polynomial is orthogonal to ψ 0 (X) = 1, which means E [ψ α (X)] = 0 ∀α = 0 and Var [ψ α (X)] = E (ψ α (X)) 2 = 1 ∀α = 0.
Thus, the time-dependent mean and standard deviation of the response can be estimated by means of a mere post-processing of the truncated PC coefficients (in Eq. ( 14)) with no additional cost as follows:

E [Y (t)] ≈ E α∈A y α (t)ψ α (X) = y 0 (t), (16)
σ 2 Y (t) = Var [Y (t)] ≈ Var α∈A y α (t)ψ α (X) = α∈A α =0 y 2 α (t). ( 17
)
3 Stochastic time-warping polynomial chaos expansions for random oscillations

Introduction

An interesting problem emerges in nonlinear oscillating systems possessing a limit cycle1 which may depend on the uncertain parameters. Limit cycle oscillations (LCO) represent a class of time-dependent problems that plays an important role in several fields, see e.g. aerospace engineering [START_REF] Bunton | Limit cycle oscillation characteristics offighter aircraft[END_REF] and mechanical engineering [START_REF] Sarrouy | Stochastic study of a non-linear self-excited system with friction[END_REF] This is achieved by minimizing the Euclidean distance between the distinct trajectories and the reference counterpart. From a similar perspective, [START_REF] Witteveen | An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase[END_REF] interpolated the oscillatory responses on the phase space to obtain in-phase oscillations. Inspired by the two mentioned approaches, a non-intrusive time transform, which consists in finding a suitable stochastic warping of the time line to increase the similarity between different trajectories in the transformed (warped) time scale, is introduced in this section. The proposed approach focuses on increasing the frequency and phase similarity of the considered trajectories in problems involving periodicity.

It is worth noting that in the engineering literature, the time-warping technique has been of interest for decades. In the context of voice recognition, [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] first proposed the time-warping to eliminate the timing differences and obtain maximum coincidences between two speech patterns. [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF] introduced a novel cost function to determine the time-warping function. Later, Ramsay and Li (1998) used the technique under the name "curve registration" for biological data. The essential idea consists in the registration (or alignment) of salient curve features by means of a suitable smooth monotone transformation of the temporal variable t. The actual analyses are then carried out on the aligned curves. Note that the same idea can also be conducted in the spatial domain. For instance, [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF] showed particular applications of registering the outcomes over surfaces or volumes in medical imaging.

Herein, we are adding one dimension to the time-warping technique by incorporating the effects of uncertainties in the transformation function. This results in a stochastic time-transform framework. Indeed, due to the inherent randomness of the stochastic problem, a time transformation function with deterministic parameters is not suitable. Therefore, stochastic transform parameters must be used and will be cast as functions of the original random parameters. The theoretical foundation of this work was originally presented by Mai and Sudret (2015b).

Stochastic time-warping polynomial chaos expansions

Consider a dynamical system (e.g. a structural dynamic or chemical system) whose behaviour is modelled by a system of ordinary differential equations (ODEs):

dy dt = f (y, ξ, t), (18)
where the initial condition is y(t = 0) = y 0 and the random vector ξ comprises independent second-order random variables defined over a probability space (Ω, F, P). ξ may include the parameters governing the system behaviour, e.g. masses, stiffness, damping ratio, reaction parameters, frequency and amplitude of excitation. The initial condition can also be uncertain, in which case it becomes a random variable belonging to ξ. The time-dependent response of the system is denoted by y(t, ξ). Without loss of generality, we consider one component of the output quantity, e.g. y(t, ξ) with the initial condition y(t = 0) = y 0 . At each time instant, y(t, ξ)

is assumed to be a second-order random variable. As in Wan andKarniadakis (2005, 2006a); [START_REF] Witteveen | An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase[END_REF]; Le Maître et al. ( 2010), herein we focus on the class of problems when y(t, ξ) is an oscillatory response with random frequencies and amplitudes.

The time-dependent response y(t, ξ) is represented by time-frozen PCEs as:

y(t, ξ) = α∈A y α (t)ψ α (ξ) + (t). (19)
A virtual time variable τ (t, ξ), which is obtained by a stochastic time-warping, is introduced as follows:

τ (t, ξ) = Nτ i=1 c i (ξ) f i (t) = F (t, ξ), (20)
where {f i (t), i = 1, . . . , N τ } are functions of time t and {c i (ξ), i = 1, . . . , N τ } are coefficients which depend on the input random variables ξ. The coefficients c i (ξ) can be represented by PCEs as:

c i (ξ) = α∈N M c iα ψ α (ξ), (21)
where ψ α (ξ) and c iα are respectively the orthonormal polynomial functions and the coefficients of the expansion. The only constraint on the time-warping is that τ is a strictly monotonically increasing function of t. Then the inverse transform may be cast as:

t(τ, ξ) = F -1 (τ, ξ). ( 22
)
Note that, in the sequel, linear transform of the form:

τ (t, ξ) = k(ξ) t + φ(ξ) (23)
are considered. Such a presentation is used to ease the presentation. For each realization ξ 0 , i.e.

each trajectory of the system response, we assume a one-to-one mapping between t and τ . The response trajectory may then be represented in the transformed (warped) time scale by:

y(τ, ξ) = β∈B y β (τ )ψ β (ξ) + (τ ), (24)
in which B is the truncation set of the multi-indices β. The inverse time transform allows one to obtain the PCEs of the response in the physical time scale as follows:

y(t, ξ) = y(F -1 (τ, ξ), ξ). ( 25
)
The objective is to find a suitable time-warping defined by Eq. ( 20) and ( 21) so that the cardinality of B remains small (i.e. low-degree PCEs can be used) to achieve an acceptable error (τ ) even at late instants. This can be obtained if the trajectories y(τ (t, ξ)) become in-phase, as suggested by Le Maître et al. ( 2010) and [START_REF] Witteveen | An alternative unsteady adaptive stochastic finite elements formulation based on interpolation at constant phase[END_REF]. First, a deterministic reference trajectory y r (t) is introduced. The stochastic time-warping (Eq. ( 20)) is determined by maximizing the similarity between y(τ (t, ξ)) and the reference counterpart y r (t) for all values of ξ, which makes the responses become in-phase. This allows the effective computation of Eq. ( 24).

Having at hand the time-warping (Eq. ( 20)) and the PCEs of the response in the virtual time line τ (Eq. ( 24)), one can finally obtain the PCEs in the physical time line t by conducting the inverse time-warping. The proposed non-intrusive time-warping approach is explained in detail in the following. For the sake of clarity, it is graphically summarized in Figure 2.

• One first chooses a reference trajectory y r (t) which is for instance obtained by considering the mean values of the input vector ξ, i.e. y r (t) = y(t, E [ξ]). In general, y r (t) may be Figure 2: Stochastic time-warping approach: computation of PCEs any realization of the response quantity y(t) obtained with a specific sample ξ 0 . For the numerical case studies considered in the current chapter, the choice of y r (t) did not affect the accuracy of the final results.

• Let us start now the time-warping, which consists in transforming the time line with the purpose of increasing the similarity between different realizations of the output y(t, ξ). Assume that one is given a set of trajectories y i (t) ≡ y(t, ξ i ), i = 1, . . . , n for n realizations of ξ corresponding to an experimental design in the input space D Ξ . Then for the realization #i, i = 1, . . . , n, the following steps are performed:

-Define a linear time-warping τ = k i t + φ i . In general, the functions f i (t) in Eq. ( 20) might be polynomials of t. However, when investigating the problem of vibration with random frequencies, a linear transform usually suffices. This is due to the periodicity of the considered response trajectories. In the intrusive time transform approach (Le Maître et al., 2010), although a linear warping function is not specified for the considered examples, the resulting transformed time τ eventually represents a linear relationship when plotted against t. [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF] also used a linear warping function. In particular, given the complexity of the problems under investigation, use of a linear function facilitates the inverse transform in the next phase, which is highly convenient. This linear warping represents two actions, namely scaling and shifting, respectively driven by the parameters k i and φ i . The time line is stretched (resp. compressed) when k i > 1 (resp. 0 < k i < 1) and is shifted to the left (resp. to the right) when φ i < 0 (resp. φ i > 0). In fact, the scaling factor k i (resp. shifting factor φ i ) allows to maximize the similarity in frequency (resp. phase) between the considered trajectories.

-Determine the parameters (k i , φ i ) governing the time-warping as the solution of an optimization problem which aims at maximizing the similarity between the response trajectory y i (k i t + φ i ) and the reference counterpart y r (t). The details of the optimization problem, in which a measure of similarity is introduced, will be described in Section 3.3.

-Represent y i (t) on the transformed time line τ . For this purpose, one chooses a grid line of τ with the desired time interval. In fact, the finer the grid is, the smaller is the error introduced by the interpolation. The trajectory y i (t) is projected onto

τ i = k i t + φ i to obtain y i (τ i ).
In order to assure that all transformed time lines τ i start at 0, when t ≤ t 0 , one uses the following transform

τ i = k i t 0 + φ i t 0 t. The small value t 0 is chosen so that k i t 0 + φ i > 0 ∀i = 1, . . . , n.
For instance, t 0 = 0.2 s is used for the numerical applications that follow. Finally the projected trajectory is linearly interpolated on the selected time line τ yielding y i (τ ).

• One builds PCEs of k(ξ), φ(ξ) and y(τ, ξ) using the realizations {k i , φ i , y i (τ )}, i = 1, . . . , n as the experimental design (or training set):

k(ξ) = γ∈G k γ ψ γ (ξ) + k , (26)
φ(ξ) = θ∈T φ θ ψ θ (ξ) + φ , (27)
y(τ, ξ) = β∈B y β (τ )ψ β (ξ) + y (τ ). ( 28
)
In the above equations, γ, θ and β are multi-indices belonging to the truncation sets G, T and B of the expansions. k γ , φ θ and y β (τ ) are coefficients computed by means of sparse adaptive PCEs [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]. k(ξ) and φ(ξ) are scalar quantities, therefore the computation of their PCE models is straightforward. However, for the vector-valued response y(τ, ξ), it might be computationally expensive when the number of discretization points of the τ -line is large. This computational cost can be reduced significantly by coupling PCEs with the principal component analysis [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF]. The combination of PCA and PCEs will be described in detail in Section 3.4.

Determination of time-warping parameters

This section describes the optimization problem used for determining the parameters k and φ of the time-warping process. We first propose a function to measure the similarity between two trajectories y 1 (t) and y 2 (t): The parameters (k i , φ i ), i = 1, . . . , n are determined as the maximizers of the similarity measure between y i (τ ) and y r (t). The objective function reads:

g(y 1 (t), y 2 (t)) = T 0 y 1 (t)y 2 (t)dt y 1 (t) y 2 (t) , (29) in which
g(k i , φ i ) = T 0 y i (k i t + φ i )y r (t)dt y i (k i t + φ i ) y r (t) . (30)
Note that the optimal warping parameters (k i , φ i ) are different for each trajectory. This results in varying total durations of the trajectories after the warping process. This also occurred in the intrusive time transform approach (Le Maître et al., 2010, Figure 4). The objective function is therefore computed on the overlapped duration between the warped trajectory and the reference one.

Let us now examine the solution (k i , φ i ) of the proposed optimization problem. The constraint that τ is a strictly monotonically increasing function of t requires that k i > 0. In case y r (t) and y i (t, ξ i ) are both monochromatic signals, the value of k i that maximizes their similarity in frequency is unique. However, there are multiple values for the shifting factor φ that make the considered trajectories in phase. This will be investigated in the next paragraph. is also a maximizer. In addition, for the sake of simplicity, it is preferable that φ is as close to 0 as possible, i.e. the time line of the scaled trajectory is shifted as least as possible. Therefore, the selected value of φ needs to satisfy the condition that the shifted distance (in time) is not larger than 1/4 of the period T r of the reference trajectory y r (t), i.e. |φ| ≤ T r /4. This constraint ensures that the solution is unique. By adopting the constraint on φ, one finds the solution (k, φ) = (2, 0) for the considered example. (

) 31
This problem can be solved by means of global optimization methods.

Principal component analysis and time-warping polynomial chaos expansions

The instant-wise application of PCEs to model the response in the transformed time line (Eq. ( 28))

might lead to an important computational burden when the discretized vector τ is of large length.

To overcome this issue, [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF] proposed a two-step approach which com-

y(τ, ξ) = ȳ(τ ) + K i=1 A i (ξ) ṽi (τ ), (32)
where ȳ(τ ) is the empirical mean vector, ṽi (τ ) is an empirical eigenvector determined with Y and A i (ξ) is a finite variance random variable. Only a few eigenvectors are retained in the decomposition, which leads to:

y(τ, ξ) = ȳ(τ ) + K i=1 A i (ξ) ṽi (τ ) + 1 (τ ). ( 33
)
The number of principal components is selected so that the relative error 1 -

K i=1 λ i K i=1 λ i
is smaller than a prescribed threshold, e.g. = 0.01. The samples of the random coefficient A i (ξ) can be obtained using a i = (Y -Ȳ ) ṽi with Ȳ = {ȳ(τ ), . . . , ȳ(τ )} being a n × K matrix obtained by replicating n times the empirical mean ȳ(τ ). The computed samples of A i (ξ) are then used as the experimental design to compute the PCE of this random coefficient:

A i (ξ) = α∈A c i,α ψ α (ξ) + 2,i . (34)
Finally, the response in the transformed time scale is represented by coupling PCA and PCEs as follows:

y(τ, ξ) = ȳ(τ ) + K i=1 α∈A c i,α ψ α (ξ) ṽi (τ ) + (τ ). ( 35
)
Note that [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF] introduced a measure of the upper bound of the total error induced by the truncation of the principal component analysis and the approximation of the random coefficients A i (ξ) by PCEs. The reader is referred to the mentioned publication for more details. Herein this error measure can be used as an indicator of the accuracy of the computed surrogate models.

Predicting random oscillations with time-warping polynomial chaos expansions

Let us now demonstrate the use of time-warping PCEs to predict responses of the model given a new set of input parameters ξ . For the sake of clarity, the procedure is depicted in Figure 4 and explained in two steps as follows:

• First, one predicts k(ξ ), φ(ξ ) and y(τ, ξ ) using the computed PCEs in equations ( 26), ( 27) and ( 35).

• Second, one maps y(τ, ξ ) into y(t, ξ ) using the inverse time-warping t = τφ(ξ ) k(ξ ) . To this end, the discretized trajectory in the warped time, say y(τ 1 , ξ ), . . . , y(τ K , ξ ) is attached to the real time instants The accuracy of the prediction #i is indicated by the relative error which reads:

t 1 = τ 1 -φ(ξ ) k(ξ ) , . . . , t K = τ K -φ(ξ ) k(ξ )
val,i = K t=1 (y(t, ξ i ) -ŷ(t, ξ i )) 2 K t=1 (y(t, ξ i ) -ȳ(t, ξ i )) 2 , ( 36
)
where ŷ(t, ξ i ) is the output trajectory predicted by PCEs and ȳ(t, ξ i ) is the mean value of the actual response time series y(t, ξ i ) which is obtained with the original numerical solver. The above formula is also used to assess the accuracy of the predicted time-dependent statistics (i.e.

mean and standard deviation).

These problems are solved in the UQLab framework [START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF], more specifically using the least angle regression algorithm implemented in the polynomial chaos expansion module [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF].

Rigid body dynamics

We first consider the rotation of a rigid body described by Euler's equations [START_REF] Peraire | 3d rigid body dynamics: equations of motion[END_REF]. The conservation of angular momentum reads:

         M x = I xx ẋ -(I yy -I zz ) y z, M y = I yy ẏ -(I zz -I xx ) z x, M z = I zz ż -(I xx -I yy ) x y, (37)
         ẋ(t) = y(t) z(t), ẏ(t) = ξ x(t) z(t), ż(t) = -x(t) y(t). (38)
The initial conditions are set equal to x(0) = 0, y(0) = 1, z(0) = 1. Assume that ξ is modelled by a random variable with uniform distribution: ξ ∼ U(-1, 1). Suppose a solver of the coupled ODEs is available. For any realization of ξ, this solver provides discretized trajectories {{x(t i ), y(t i ), z(t i )} , t i = 0, ∆ t , . . . , K ∆ t ≡ T }. In this example, the equations are solved using the Matlab ordinary differential equation solver ode45 (Runge-Kutta method, total duration T = 50 s, time step ∆ t = 0.01). We aim at building PCEs of the angular velocity x(t) as a function of the random variable ξ. Note that the corresponding polynomial functions are from the family of orthonormal Legendre polynomials since ξ is uniformly distributed.

Figure 5 depicts a set of 50 trajectories of x(t) obtained for different realizations of the random variable ξ. This set is used as the experimental design for fitting the time-frozen PCEs.

x(t) are oscillatory trajectories which fluctuate around zero at different frequencies. This is a typical example of the problem of stochastic oscillation with uncertain frequencies (Wan andKarniadakis, 2005, 2006a). At the early instants (t < 10 s), one can differentiate between the distinct trajectories, whereas this is hardly the case at later instants, since the patterns are mixed up completely. Due to the growing difference in frequency and phase, x(t, ξ) is more and more non-linear as a function of ξ for increasing t (Figure 6(A)). Subsequently, the probability density function of x(t) becomes bi-modal at late instants (Figure 6(B)). This explains why increasingdegree time-frozen PCEs are required in order to represent x(t) properly. As analyzed previously, this is not a sustainable approach since the required degree of PCEs will certainly become too high at some point. Time-frozen sparse PCEs are now utilized to model the variability of the response trajectories, and exemplify the deficiency of such an approach. At each instant t, an adaptive PCE scheme with candidate polynomials up to total degree 20 is used (Eq. ( 19)) based on the available 50 data points from the experimental design made of the 50 trajectories. The PCE model which results in the smallest leave-one-out (LOO) error is retained. Figure 7 depicts the LOO error of these time-frozen PCEs, which is increasing in time, showing that the accuracy of the PCE model degenerates. For validation purpose, a set of 10, 000 trajectories is computed using the ode45 Matlab solver. Figure 8 depicts two particular response trajectories predicted by time-frozen PCEs versus the actual responses obtained by numerically solving the system of ordinary differential equations (38). After 15 s (when the LOO error is approximately 10 -2 ) the PCE prediction deviates significantly from the actual trajectory. In particular, there are signs of instability in the PCE model, e.g. the PCE-based prediction for consecutive instants differ noticeably in terms of accuracy.

We now consider the time-dependent mean and standard deviation of the response x(t) which are depicted in Figure 9. In the early time instants (t < 15 s), time-frozen PCEs represent the statistics with relatively small error compared to Monte Carlo simulation (MCS). However, after 15 s, the accuracy declines quickly. In particular, PCEs cannot mimic the oscillatory behaviour of the standard deviation. Another interpretation is that even degree-20 time-frozen PCEs cannot capture the complex distribution of the response at late time instants.

Let us now apply the time-warping approach to pre-process the trajectories x(t). Provided that the initial condition is equal to 0, it suffices to use a linear time-warping τ = k t. For each computed realization of the angular velocity x(t, ξ i ), i = 1, . . . , 50, the parameters k i is estimated as the maximizer of the similarity measure described in Eq. ( 30). Note that the same 50 trajectories are used as the experimental design for this approach and the reference trajectory The relative LOO error is 3.82 × 10 -4 , which indicates a high accuracy of the PCE model.

The time-warping is carried out using the estimated parameters and the responses are interpolated into the transformed time line τ , leading to in-phase trajectories x(τ ) (see Figure 10(A)).

As expected, x(τ ) are smooth functions of ξ at all instants, which allows the effective use of PCEs (Figure 10(B)).

Principal component analysis (PCA) is then conducted on the obtained transformed trajec- Using PCE to represent the first PCA coefficient, the obtained relative LOO error is 7.7 × 10 -3 .

It is also clear that it is harder to represent the higher mode PCA coefficients by PCEs, as was observed by [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF]. However, it is worth noting that most of the stochastic features of the response is captured by the first few components.

Kraichnan-Orszag model

Let us investigate dynamical systems with random initial conditions, e.g. the so-called Kraichnan-Orszag three-mode problem. It was introduced by [START_REF] Kraichnan | Direct-interaction approximation for a system of several interacting simple shear waves[END_REF] to model a system of several interacting shear waves and later was studied by [START_REF] Orszag | Dynamical properties of truncated Wiener-Hermite expansions[END_REF] in the case of Gaussian initial conditions. This model is described by the following system of ODEs:

         ẋ(t) = y(t) z(t), ẏ(t) = z(t) x(t), ż(t) = -2 x(t) y(t). ( 39
)
The initial condition of x(t) is considered stochastic, i.e. x(t = 0) = α + 0.01 ξ with ξ ∼ U[-1, 1] whereas y(t = 0) = 1.0, z(t = 0) = 1.0. Herein, we consider α = 0.99 as investigated by [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] with the time-dependent PCEs. Note that when α is in the range [0, 0.9], the responses are insensitive to the initial conditions. For α ∈ [0.9, 1], there is a strong dependence of the responses on the initial state. Figure 13(A) depicts the large discrepancies between time-histories of x(t) due to a minor variability of the initial condition x(t = 0). The surrogate model of the response x(t) is computed with time-frozen and time-warping PCEs using an experimental design of size N = 50 (Figure 13(A)). On the one hand, adaptive sparse PCEs with candidate bases up to total degree 20 are used for the time-frozen approach.

On the other hand, a time-transform scheme τ = k t with one governing parameter is used for the time-warping scheme. The trajectories resulting from the time-warping process are depicted in Figure 13(B). The adaptive sparse PCE representing k has the relative LOO error 2.2×10 -6 . The first 13 principal components are retained so that 99.9% of the response's variance is explained.

The relative LOO errors of PCEs for the first two components are 9.4 × 10 -5 and 7 × 10 -3 , respectively.

The time-warping PCE model is then validated by accessing the accuracy of its predictions.

Figure 14 plots two specific predictions of the surrogate model which are graphically indistin-guishable from the actual time-histories obtained with the original Matlab solver. Only 1.27% of the 10, 000 predictions experiences a relative error larger than 0.1. Regarding the mean and standard deviation trajectories (Figure 15), the time-warping approach leads to respective relative errors 2.1×10 -4 and 5.3×10 -4 , which shows an excellent agreement between the predictions and the references. These figures also show that the time-frozen sparse PCEs computed with the same experimental design of size 50 lead to predictions which are not sufficiently accurate. This numerical application illustrates the potential application of the proposed time-warping approach to systems subject to uncertain initial conditions. The excellent performance of the approach is even more impressive given the chaotic behaviour of the considered system, i.e. the responses are strongly sensitive with respect to a minor variability of the initial condition.

Duffing oscillator

Let us consider a non-linear damped single-degree-of-freedom (SDOF) Duffing oscillator under free vibration, which is described by the following equation of motion:

ÿ(t) + 2 ω ζ ẏ(t) + ω 2 (y(t) + y 3 (t)) = 0. ( 40
)
The oscillator is driven by uncertain parameters ξ = (ζ, ω, ) described in Table 1. The initial conditions are considered deterministic with y(t = 0) = 1 and ẏ(t = 0) = 0. Note that a simplified form of this equation which represents an undamped linear oscillator was used in other publications for illustrating the time-dependent generalized polynomial chaos [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF], the intrusive time-transform approach (Le Maître et al., 2010) and the flow map composition PCEs [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF]. relative LOO errors 1.87 × 10 -5 and 2.08 × 10 -4 respectively, which indicates a high level of accuracy. PCA is then applied to retrieve eight principal components that results in the PCA truncation error smaller than 1 × 10 -3 . The relative LOO errors of PCE models for the first two components are 8 × 10 -4 and 4 × 10 -3 , respectively.

An independent validation set of 10, 000 runs is used to judge the accuracy of the PCE models. In terms of time-dependent statistics (Figure 18), time-frozen PCEs can predict rather well the mean trajectory, however fail to represent the standard deviation after early instants (t > 3 s).

In contrast, the time-warping approach provides excellent accuracy on the mean and standard deviation time histories. The relative discrepancies between mean and standard deviation time histories predicted by time-warping PCEs with the reference trajectories are 3.27 × 10 -5 and 3.47 × 10 -4 , respectively.

Oregonator model

We consider now the Oregonator model which describes the dynamics of a well-stirred, homogeneous chemical system governed by a three species coupled mechanism. Note that this benchmark problem was used by Le Maître et al. (2010) to illustrate the intrusive time-transform approach.

This chemical system undergoes an oscillation governed by the following system of ODEs:

         ẋ(t) = k 1 y(t) -k 2 x(t) y(t) + k 3 x(t) -k 4 x(t) 2 , ẏ(t) = -k 1 y(t) -k 2 x(t) y(t) + k 5 z(t), ż(t) = k 3 x(t) -k 5 z(t), (41)
in which (x, y, z) denotes the three species concentration and the coefficients k i , i = 1, . . . , 5 are the reaction parameters. Hereafter, all the reaction parameters are considered independent random variables with uniform and normal distributions (see Table 2). It is worth noting that Le Maître et al. ( 2010) considered only k 4 and k 5 as uniform random variables while fixing the remaining parameters (i.e. k 1 = 2, k 2 = 0.1, k 3 = 104). The initial condition is (x 0 , y 0 , z 0 ) = (6, 000; 6, 000; 6, 000), which corresponds to a deterministic mixture. We aim at building PCEs of the concentration x(t) as a function of the random parameters ξ = (k 1 , k 2 , k 3 , k 4 , k 5 ). We now apply the proposed non-intrusive time-warping approach to this problem. Note that only 50 trajectories of x(t) are used as an experimental design for this approach. A linear time- Finally, the time-warping PCE scheme is applied to surrogate the responses y(t) and z(t) of the system using the same experimental design of size 50 and the same procedure. Figure 22 shows a great agreement between two specific trajectories, the mean and standard deviation of (x, y, z) in the state-space predicted by time-warping PCEs and the reference functions.

Forced vibration of a Bouc-Wen oscillator

In the previous case studies, self-oscillating systems were considered. In this example, we show that the proposed approach is also applicable to forced-vibration systems. Let us now consider the SDOF Bouc-Wen oscillator [START_REF] Kafali | Seismic fragility analysis: Application to simple linear and nonlinear systems[END_REF] subject to a stochastic excitation.

The equation of motion of the oscillator reads: in which ζ is the damping ratio, ω is the fundamental frequency, ρ is the post-to pre-yield stiffness ratio, γ, α, β, n are parameters governing the hysteretic loops and the excitation x(t) is a sinusoidal function given by x(t) = A sin(ω x t).

   ÿ(t) + 2 ζ ω ẏ(t) + ω 2 (ρ y(t) + (1 -ρ) z(t)) = -x(t) , ż(t) = γ ẏ(t) -α | ẏ(t)| |z(t)| n-1 z(t) -β ẏ(t) |z(t)| n . ( 42
Deterministic values are used for the following parameters of the Bouc-Wen model: ρ = 0, γ = 1, n = 1, β = 0. The remaining parameters ξ = (ζ, ω, α, A, ω x ) are considered independent random variables with associated distributions given in Table 3.

One aims at representing the oscillator displacement y(t) as a function of the uncertain input parameters using time-frozen and time-warping PCEs. To this end, 100 simulations of the oscillator are carried out using the Matlab solver ode45 with time increment ∆ t = 0.005 s for the total duration T = 30 s and initial condition y(t = 0) = 0, ẏ(t = 0) = 0. The displacement trajectories are depicted in Figure 23(A). First, the time-frozen sparse PCEs are computed with candidate polynomials up to total degree 20. For this case study, a time-warping scheme τ = k t with only one parameter is used.

After the time-warping process, the trajectories become in-phase as depicted in Figure 23(B).

Adaptive sparse PCE representing k has the relative LOO error 5 × 10 -5 . In order to achieve a truncation error 1 smaller than 1 × 10 -3 , 13 first principal components are retained in PCA.

The relative LOO errors of PCEs for the first two components are 6 × 10 -3 and 6.21 × 10 -2 , respectively.

Let us validate the accuracy of the time-warping PCE model. In Figure 24, two specific predictions of the PCE model are plotted against the actual responses obtained with the original Matlab solver. A remarkable agreement can be observed. Among 10, 000 validations, only 4.87% has a relative error larger than 0.1. Regarding the time-dependent mean and standard deviation of the oscillator, time-warping PCE-based estimates outstandingly match the reference trajectories (Figure 25). Only a minor discrepancy can be observed at the end of the considered time duration T = 30 s, which is due to the modest number of simulations used as the experimental design. The corresponding relative errors are both 2.4 × 10 -3 . On the contrary, time-frozen PCEs exhibit a low level of accuracy after 5 seconds. It is worth noting that in the current case study, we considered both uncertainties from the mechanical properties and the excitations. In particular, complicated hysteretic behaviour was investigated. To the best of the authors' knowledge, this is the first time that such a system is considered in the literature of uncertainty quantification for the purpose of deriving timedependent surrogate models.

Discussion

The various numerical applications in chemical and mechanical engineering have proved the effectiveness of the time-warping PCE approach, which may be shortly explained as follows. It was observed that when represented in the space of the temporal variable t, the system's responses are increasingly non-linear functions of the uncertain parameters. When projecting the responses onto a suitable space, in this case the transformed time line τ , the resulting trajectories become smooth functions of the uncertain input parameters, whose complexity does hardly increase with time. Therefore, PCEs can be applied effectively to the projected responses and represent well the solutions at late instants. In this paper, a measure of similarity was proposed to define a suitable space for projecting the responses, which exploits the periodicity of the trajectories.

Further investigations are required to clearly determine such a suitable space in a more general case.

In the proposed approach, the virtual time τ is a function of the uncertain parameters ξ. In other words, the basis "function" τ onto which the responses are projected is not deterministic. This differs significantly from approaches commonly used in the literature, in which the response trajectories are first projected onto a set of deterministic reduced basis determined a priori using a set of numerical simulations of the system. This is usually done with a simple linear transform, for instance data compression techniques such as principal component analysis or wavelet decomposition.

When analysing further, one discovers a particular feature which constitutes a major difference between the classical time-frozen PCE approach and the proposed time-warping method.

The PC coefficients y β (τ ) in the time-warping representation (Eq. ( 24)) are functions of τ , therefore being dependent on ξ. This contradicts the representation of time-frozen PCEs (Eq. ( 19)), in which t and ξ intervene in the solution in a separated manner.

From a more general perspective, the effectiveness of the approach can be explained by analysing the functionalities of the time-warping process and PCEs. The most important feature of an oscillatory trajectory consists in its spectral content, which is characterized by the vibration periodicity. The other feature is the temporal content characterized by the vibration amplitude. The pre-processing step handles partially the dynamics of the system by dealing with the frequency content. Using the time-warping process, the resulting trajectories have similar frequencies and phases. In other words, in terms of frequencies, the transformed trajectories exhibit a similar dynamical behaviour, which is close to that of the reference trajectory. The other aspect of the dynamics, i.e. the random temporal amplitude of the trajectories, is handled with sparse PCEs. As a summary, the dynamics is captured by the time-warping process, whereas the uncertainties are represented by PCEs.

As explained, sparse PCEs alone are not capable of dealing with the dynamics. The proposed approach illustrates a novel way to solve stochastic dynamical problems, in which a specialized technique might be used to capture the dynamical aspect whereas sparse PCEs are used to propagate uncertainties. From this perspective, [START_REF] Yaghoubi | Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation[END_REF] have recently applied the warping-based approach in the frequency domain to surrogate the frequency response function of mechanical systems. This principle is further developed by [START_REF] Mai | Surrogate modelling for stochastic dynamical systems by combining NARX models and polynomial chaos expansions[END_REF] to tackle more complex problems in which non-linear uncertain structures subject to stochastic motions are of interest and where the response trajectories are non-stationary, i.e. they do not show pseudoperiodic oscillations. The projection of the responses onto a special basis made of auto-regressive functions will allow us to represent the non-linear dynamical behaviour of the systems.

Finally, it is worthwhile mentioning that the proposed methodology is fully non-intrusive, i.e. the surrogate models of the systems' response trajectories are obtained by using a pre-computed set of trajectories related to an experimental design. In this respect, the methodology is readily applicable to any other problems featuring randomized limit cycle oscillations.

Conclusions and perspectives

Polynomial chaos expansions (PCEs) represent an effective metamodelling technique which has been efficiently used in several practical problems in a wide variety of domains. It is, however, well known that PCEs fail when modelling the stochastic responses at late instants of dynamical systems. In this paper, we pointed out the cause of the failure, which is mainly associated with the large dissimilarities between distinct responses introduced by the variability of the uncertain parameters.

To address the above issue, we suggested an approach which consists in representing the responses into a virtual time line where the similarities between different response trajectories are maximized. The virtual time line is obtained by warping, i.e. scaling and shifting, the original time grid. The parameters governing the trajectory-dependent time warping are determined by means of a global optimization problem using an objective function herein introduced to quantify the similarity between distinct trajectories. The proposed approach allows one to effectively solve complex benchmark problems from mechanics and chemistry using only low-order PCEs. This approach also suggests that when representing the original response quantities onto a suitable transformed space, the complexity of the responses may reduce significantly, thus allowing more effective application of PCEs. In general, pre-processing the experimental design before applying PCEs is a promising approach that needs further investigation.

  represented the dynamic response trajectories as functions of the phase φ instead of time t in order to obtain in-phase vibrations. The phases are extracted from the observations, based on the local extrema of the time series. The response trajectories are then transformed from time-histories to phase-histories. PCEs are eventually applied in the phase space. Finally, Le Maître et al. (2010) represented the responses in a rescaled time τ such that the dynamic responses vary in a small neighbourhood of a reference trajectory. The time scale τ is intrusively adjusted at each time step so that the distance between the dynamic response and the reference solution is minimized, thus in-phase vibrations are achieved.

  among others. Use of PCEs to represent LCO systems has attracted a large attention and actually almost all novel ideas with PCEs are applied first to LCO systems or systems involving periodicity. For instance, Wan and Karniadakis (2006a) used multi-element PCEs whereas Beran et al. (2006b) proposed different methods namely use of Haar wavelets as local bases or use of B-spline functions. These approaches aim at resolving the highly nonlinear behaviour of LCO responses in the stochastic domain. There are also techniques that are designed specifically for LCO. Le Maître et al. (2010) proposed an intrusive time transform of the trajectories which aims at representing the transformed time-histories in a small neighbourhood of a reference trajectory, i.e. to reduce their variability by making them in-phase. A transformed time line τ is introduced, of which the varying clock speed τ = dτ dt is adjusted in an intrusive setup at each time step.

  t)y 2 (t)dt is the inner product of the two considered time histories and • is the associated L 2 -norm. In practice, the trajectories are discretized and thus, the inner product (resp. the L 2 -norm) becomes the classical dot product between two vectors (resp. the Euclidean norm). By the Cauchy-Schwarz inequality, this similarity measure always takes values in the interval [0, 1]. It attains its maximum when the considered trajectories have the same frequency and phase content.

Figure 3

 3 Figure3depicts the objective function g(k, φ) as a similarity measure between the reference trajectory y r (t) = sin(π t) and a response y(t) = sin(2 π t). The two trajectories are chosen in such a way that (k, φ) = (2, 0) is the maximizer of g(k, φ). However, there are three global maxima in the depicted interval[-1.5, 1.5] of φ. This is due to the fact that in the virtual time line τ , if the transformed trajectory y(τ ) is shifted (whether to the left or to the right) a distance equal to one half of the period T r = 2 s of the reference counterpart, the similarity measure reaches another global maximum. In fact, if T r /4 ≤ φ ≤ T r /2 (resp. -T r /2 ≤ φ ≤ -T r /4) maximizes the similarity measure, then φ -T r /2 (resp. φ + T r /2) in the interval [-T r /4, T r /4]

Figure 3 :

 3 Figure 3: Similarity measure as a function of k and φ

  bines principal component analysis (PCA) and PCEs. The first step consists in conducting PCA to capture the stochastic features of the random vector-valued response with a small number of deterministic principal components and the associated non-physical random variables. The second step relies on representing the resulting random variables with adaptive sparse PCEs.Consider a sample set of the response trajectories Y = y (1) (τ ), . . . , y (n) (τ ) represented at the discretized points {τ 1 , . . . , τ k } in the transformed time line. By stacking up the discretized responses, one obtain a matrix of trajectories of size n × K denoted by Y . The response can be represented by PCA as follows:

Figure 4 :

 4 Figure 4: Stochastic time-warping approach: prediction of the response trajectories using PCEs

  in which M x , M y , M z are the external moments, I xx , I yy , I zz are the moments of inertia and x, y, z are the angular velocities about the principal axes. In the case when the rigid body rotates freely under no external excitation, i.e. M x = M y = M z = 0 and I xx = 1ξ 2 I yy , I zz = 1 + ξ 2 I yy , one obtains the following set of reduced equations:

Figure 5 :

 5 Figure 5: Rigid body dynamics -N = 50 different trajectories x(t) in the original time scale t.

  Probability density function of x(to) at to = 5, 10, 15, 20 s

Figure 6 :

 6 Figure 6: Rigid body dynamics -x(t, ξ) as a function of ξ for particular instants and its probability density function.

Figure 7 :

 7 Figure 7: Rigid body dynamics -Leave-one-out error of time-frozen PCEs.

Figure 8 :

 8 Figure 8: Rigid body dynamics -Two particular trajectories and their predictions by time-frozen and time-warping PCEs.

Figure 9 :

 9 Figure 9: Rigid body dynamics -Mean and standard deviation of the trajectories: comparison of the two approaches.

  Figure 10: Rigid body dynamics -Different trajectories x(τ ) in the warped time scale τ and x(τ ) as a function of the random variable ξ.

Figure 11 :

 11 Figure 11: Rigid body dynamics -The first eight principal components.

Figure 12 :

 12 Figure 12: Rigid body dynamics -PCA truncation-induced error 1 and PCE approximation error 2 normalized by trace( Σ) with Σ being the empirical covariance matrix of the sample set Y and the upper bound = √ 1 + √ 2 2 of the total error.

Figure 8

 8 Figure8depicts two specific realizations of the angular velocity x(t) predicted by timewarping PCEs, which are plotted together with the predictions by time-frozen sparse PCEs and the actual responses obtained by the numerical solver. As mentioned previously, one observes that starting from 15 s, the direct approach encounters instability, which results in inaccurate predictions. The time-warping approach allows one to improve notably the quality of the surrogate model. The predictions by time-warping PCEs are in excellent agreement with the actual responses. A relative error exceeding 0.1 is recorded in only 79 simulations among 10, 000 validations.In Figure9, the time-dependent mean and standard deviation of the response are plotted. Time-frozen PCEs allow one to represent the mean trajectory with relatively small discrepancy compared to the trajectory obtained with the MCS. It can faithfully predict the standard deviation at the early instants t < 15 s, however becomes suddenly unstable afterwards. In contrast, time-warping PCEs provide estimates of the statistics that are almost indistinguishable from the MCS estimates. The relative errors between the reference and predicted mean and standard deviation are 7.31 × 10 -4 and 7.19 × 10 -4 , respectively.

  Figure 13: Kraichnan-Orszag model -N = 50 different trajectories in the original and warped time scales.

Figure 14 :

 14 Figure 14: Kraichnan-Orszag model -Two particular trajectories and their predictions by timewarping PCEs.

Figure 15 :

 15 Figure 15: Kraichnan-Orszag model -Mean and standard deviation of the trajectories: comparison of the two approaches.

  Figure 16(B). A linear time-warping with two parameters, i.e. τ = k t + φ, is used for each trajectory. Using sparse PCEs of degree up to 20, the metamodels of k and φ are obtained with

Figure 17 Figure 16 :

 1716 Figure17presents two specific realizations of the displacement y(t) obtained with two distinct sets of parameters (ζ, ω, ). Without time-warping, PCEs are capable of predicting the response at the early time instants (t < 3 s), then their accuracies degenerate with time, resulting in incorrect predictions. By introducing the time-warping of the trajectories, PCEs can faithfully capture the damped oscillatory behaviour. Only 0.18% of 10, 000 predictions exhibits a relative

Figure 17 :

 17 Figure 17: Duffing oscillator -Two particular trajectories and their predictions by time-frozen and time-warping PCEs.

Figure 18 :

 18 Figure 18: Duffing oscillator -Mean and standard deviation of the trajectories: comparison of the two approaches.

Figure 19 (

 19 Figure19(A) depicts 50 trajectories among 500 realizations of x(t), which are used as the experimental design for fitting time-frozen PCEs. One notices that after 5 seconds, the different trajectories are completely out-of-phase. Time-frozen sparse PCEs with candidate polynomials up to total degree 20 are used. The PCE model actually starts degenerating at t = 3 s. In

  Figure 19: Oregonator model -N = 50 different trajectories of the response x(t). The figures in the first row are zoomed in the range [15, 20] to obtain the figures in the second row.

  transform τ = k t + φ is again utilized. The parameters k and φ are determined and sparse PCEs of k and φ are then computed. The relative LOO errors of the PCE models for k and φ are respectively 4.42×10 -5 and 4.8×10 -2 , which indicate a high accuracy. The response trajectories are interpolated into the transformed time line τ (Figure 19(B)) and adaptive sparse PCEs with candidate polynomials up to total degree 20 combined with PCA are then used. The first 18components are retained in PCA to obtain a truncation error 1 smaller than 1 × 10 -2 . The PCEs for the first two coefficients have relative errors 7.57 × 10 -4 and 1.5 × 10 -3 , respectively.A validation set of 10, 000 trajectories is used to get reference trajectories of the concentration x(t).

  Figure 20 depicts two particular realizations computed by the numerical solver (Matlab ordinary differential equation solver ode45, using a time step ∆ t = 0.01 for the total duration T = 40 s) and predictions by PCEs with and without time-warping. It is shown that without time-warping, PCEs fail to capture the oscillatory behaviour of the response. In contrast, the use of time-warping allows PCEs to predict the response with great accuracy. Only 1.24% of the predictions (among 10, 000 samples) has a relative error larger than 0.1. Figure 21 depicts the statistics of x(t) predicted by time-frozen and time-warping PCEs in comparison with MCSbased trajectories. Without time-warping, the estimates by PCEs differ significantly from the reference trajectories already from 3 s. The discrepancies then quickly increase in time. For instance, PCEs without time-warping estimate a decreasing trend in time for the standard deviation, whereas the latter actually oscillates around a constant value (around 1400) with high frequency. By introducing the time-warping pre-processing, one can use sparse PCEs to capture the complex behaviour of the time-dependent statistics of the response all along the trajectories. The relative error for the mean and standard deviation trajectories are 3.11 × 10 -4 and 3.6 × 10 -3 , respectively.

Figure 20 :

 20 Figure20: Oregonator model -Two particular trajectories x(t) and their predictions by time-frozen and time-warping PCEs. The figures in the first row are zoomed in the range[15, 20] to obtain the figures in the second row.

Figure 21 :

 21 Figure 21: Oregonator model -Mean and standard deviation of x(t): comparison of the two approaches.

Figure 22 :

 22 Figure22: Oregonator model -Trajectories of (x(t), y(t), z(t)) predicted by time-warping PCEs vs. the reference trajectories.

Figure 23 :

 23 Figure 23: Bouc-Wen oscillator -N = 100 different trajectories of the solution in the original time scale t and in the transformed time line τ .

Figure 24 :

 24 Figure24: Bouc-Wen oscillator -Two particular trajectories and their predictions by the two approaches.

Figure 25 :

 25 Figure 25: Bouc-Wen oscillator -Mean and standard deviation of the trajectories: comparison of the two approaches.

Table 1 :

 1 Uncertain parameters of the Duffing oscillator

	Parameters Distribution Mean Standard deviation Coefficient of variation ζ Uniform 0.03 0.015/ √ 3 0.2887 ω Uniform 2 π π/ √ 3 0.2887 Uniform -0.5 0.25/ √ 3 0.2887

Hereafter, we aim at building PCEs

Table 2 :

 2 Reaction parameters of the Oregonator model

	Parameters Distribution Mean Standard deviation Coefficient of variation k 1 Uniform 2 0.2/ √ 3 0.0577 k 2 Uniform 0.1 0.005/ √ 3 0.0289
	k 3 k 4 k 5	Gaussian Uniform Uniform	104 0.008 26	1.04 4 × 10 -4 / √ 2.6/ √ 3	3	0.01 0.0289 0.0577
	particular, Figure 20 shows that when used for predicting the responses, time-frozen PCE provide
	negative values of the concentration at some instants, which is non physical for the considered
	problem.

  )

						Reference							Reference
		8000				Time-frozen PCE			8000				Time-frozen PCE
						Time-warping PCE						Time-warping PCE
		6000								6000
	x(t)	4000							x(t)	4000
		2000								2000
		0								0
		0	10		20 t (s)	30		40		0	10		20 t (s)	30	40
	(A) ξ = (1.8970, 0.1001, 104.3676, 0.0077, 25.5417)	(B) ξ = (1.9481, 0.0999, 102.7929, 0.008, 27.6482)
						Reference							Reference
		8000				Time-frozen PCE			8000				Time-frozen PCE
						Time-warping PCE						Time-warping PCE
		6000								6000
	x(t)	4000							x(t)	4000
		2000								2000
		0								0
		15	16	17	t (s)	18	19	20		15	16	17	t (s)	18	19	20
	(C) ξ = (1.8970, 0.1001, 104.3676, 0.0077, 25.5417)	(D) ξ = (1.9481, 0.0999, 102.7929, 0.008, 27.6482)
	(zoom)							(zoom)

Table 3 :

 3 Uncertain parameters of the Bouc-Wen model

	Parameters Distribution Mean Standard deviation Coefficient of variation
	ζ			Uniform		0.02	0.002	0.1
	ω			Uniform		2 π	0.2 π	0.1
	α			Uniform		50	5	0.1
	A			Uniform		1	0.1	0.1
	ω x			Uniform		π	0.1 π	0.1
	0.1
	0
	-0.1
	y(t)
	-0.2
	-0.3
	0 -0.4	5	10	15 t (s)	20	25	30

Limit cycle is a closed isolated trajectory in the phase-space of self-oscillated oscillators. The nearby trajectories can either spiral in toward or away from the limit cycle.

The term "time-frozen sparse adaptive PCEs" refers to the instantaneous computation of sparse adaptive PCEs.