Next Article in Journal
A Generalized Alternating Linearization Bundle Method for Structured Convex Optimization with Inexact First-Order Oracles
Next Article in Special Issue
p-Refined Multilevel Quasi-Monte Carlo for Galerkin Finite Element Methods with Applications in Civil Engineering
Previous Article in Journal
Success History-Based Position Adaptation in Fuzzy-Controlled Ensemble of Biology-Inspired Algorithms
Previous Article in Special Issue
Uncertainty Propagation through a Point Model for Steady-State Two-Phase Pipe Flow
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Application of Generalized Polynomial Chaos for Quantification of Uncertainties of Time Averages and Their Sensitivities in Chaotic Systems

by
Kyriakos Dimitrios Kantarakias
* and
George Papadakis
Department of Aeronautics, Imperial College London, London SW7 2AZ, UK
*
Author to whom correspondence should be addressed.
Algorithms 2020, 13(4), 90; https://doi.org/10.3390/a13040090
Submission received: 3 February 2020 / Revised: 25 March 2020 / Accepted: 8 April 2020 / Published: 13 April 2020

Abstract

:
In this paper, we consider the effect of stochastic uncertainties on non-linear systems with chaotic behavior. More specifically, we quantify the effect of parametric uncertainties to time-averaged quantities and their sensitivities. Sampling methods for Uncertainty Quantification (UQ), such as the Monte–Carlo (MC), are very costly, while traditional methods for sensitivity analysis, such as the adjoint, fail in chaotic systems. In this work, we employ the non-intrusive generalized Polynomial Chaos (gPC) for UQ, coupled with the Multiple-Shooting Shadowing (MSS) algorithm for sensitivity analysis of chaotic systems. It is shown that the gPC, coupled with MSS, is an appropriate method for conducting UQ in chaotic systems and produces results that match well with those from MC and Finite-Differences (FD).

1. Introduction

Over the past few decades, modern computational methods have been very successful in predicting the evolution of systems of great complexity. However, in real–life applications, the system behavior is significantly affected by parametric uncertainties. As an example, geometric uncertainties in the shape of the front wing, diffuser, and the tip of a Formula (1) can affect its performance [1]. A small selection of other examples can be found in References [2,3,4,5].
Various techniques for the quantification of uncertainties have been developed. These include sampling techniques, such as the Monte-Carlo [6], Latin Hypercube sampling [7], and quasi Monte–Carlo [8]. However, despite their accuracy and ease of application, these methods become computationally intractable for large, real world problems, such us high fidelity simulations of turbulent flows. Other approaches, based on the spectral representation of uncertain quantities, have also been developed. These techniques are referred to as Polynomial Chaos Expansion (PCE) methods. They were first developed in Reference [9] for Gaussian uncertainties, and, later in Reference [10], they were generalized for all Probability Density Functions (PDFs) in the Askey scheme. The ability of generalized PCE (gPC) in quantifying uncertainties has been amply demonstrated in many applications, for example, in Computational Fluid Dynamic (CFD) systems [11,12,13]. Recent developments in Uncertainty Quantification (UQ) include using the PCE in data-driven approaches [14]. Spectral approaches are computationally efficient for a small number of uncertain parameters, m, but their cost increases exponentially with m (curse of dimensionality).
There are two major approaches to applying gPC, the intrusive (iPC) [5,15] and the non-intrusive (niPC) [16]. In the intrusive approach, the state vector is written in spectral form using the gPC, and a new set of state equations, the iPC equations, are extracted using Galerkin projection. On the other hand, the non-intrusive approach is based on the spectral representation of a Quantity of Interest (QoI), such as the lift of an airfoil, and the spectral coefficients are calculated by repeated evaluations of the QoI at specific quadrature locations.
For general unsteady systems, the iPC approach is not feasible for conducting UQ because the accuracy of the spectral representation of the state variables starts to deteriorate after the system starts to evolve from the initial condition. For this reason, we used the niPC approach and focused on uncertainties of time-averaged quantities and their sensitivities with respect to input parameters. Time-averaged quantities are of interest to many engineering applications, for example, in aeronautics (lift and drag around an airfoil), energy generation (fuel consumption and pollutant formation), process industries (average heat transfer, mixing rate), etc.
Although time-averages of state variables can be easily computed by collecting statistics over long time integration, the evaluation of their sensitivities to input parameters is much more involved. The adjoint sensitivity analysis approach [17] is based on a linearization of the evolution equation with respect to one or more system parameters. This approach, however, breaks down in chaotic systems, where the state variables display extreme sensitivity to minute variations in input conditions [18]. In such cases, the adjoint variables diverge; this is popularly known as the ‘butterfly effect’, and it is due to the presence of one or positive Lyapunov exponents.
The Least-Squares Shadowing (LSS) approach was developed in Reference [19] and resolves this problem. It is based on the shadowing lemma for ergodic and uniformly hyperbolic dynamical systems [20]. Although computationally expensive, the method can give values for the sensitivities of time-averaged quantities that match well with finite differences. In this work, the preconditioned Multiple-Shooting (MSS) [21,22], a computationally more efficient version of the LSS, is used to evaluate sensitivities of time-averages.
In this paper, it is shown that the non-intrusive approach to the gPC is suitable for UQ of time-averages of chaotic systems and their sensitivities, with respect to system parameters. Information on the latter, i.e., UQ of sensitivities, is important in robust design applications, where the expectation of a QoI is minimized under uncertainty. The present work follows on from previous work of the authors in this area [23].

2. Sensitivity Analysis of Chaotic Systems Using MSS

In this work, we consider a non-linear dynamical system governed by the equation
d u d t = f ( u , s ) ,
where u denotes the state vector, and s a set of design parameters. This is a set of ordinary differential equations (ODEs) that arise after spatial discretisation of partial differential equations (PDEs). The QoI is the time average of some function G ( u , s ) of the state variables and design parameters, as in
G ¯ ( s ) = lim T 1 T 0 T G ( u , s ) d t .
An example QoI from the field of aerodynamics is the time-average of the lift or the drag coefficient of an airfoil.
The conventional adjoint sensitivity analysis to determine d G ¯ ( s ) d s requires the formation of the augmented Lagrangian
Λ = 1 T 0 T G ( u , s ) + λ T d u d t f ( u , s ) d t + μ T u ( 0 ) u 0 .
After differentiation with respect to s and integration by parts, we arrive at the adjoint system
d λ d t = f u T λ + G ( u , s ) s ,
which is integrated backwards in time with terminal condition λ ( T ) = 0 . As mentioned in the Introduction, if the system in Equation (1) is chaotic and has at least one positive Lyapunov exponent, then marching Equation (4) backwards will result in an exponentially growing adjoint state vector λ .
To resolve the problem of diverging trajectories for s and s + δ s , we formulate a minimization problem that keeps the two trajectories close in phase space. To achieve this, we relax the requirement of fixed initial condition for the two nearby values of s . However, for ergodic systems, this change does not affect the time-average values and their sensitivities. Below, we state the minimization problem for a single parameter s, but an adjoint version for multiple parameters can be also formulated.
min v ( t i + ) i = 0 K v ( t i + ) 2 2 ,
subject to v ( t i + ) = v ( t i ) ,
d v d t f u v f s η f = 0 ,
f ( u ( t ) , s ) , v ( t ) = 0 ,
where t i < t < t i + 1 ( i = 0 , . . . , K 1 ) ,
where
v ( t ) = lim δ s 0 u ( τ ( t ) ; s + d s ) u ( t ; s ) δ s η ( t ) = d d s d τ ( t ) d t 1 .
The objective to be minimized, given in Equation (5a), is the L 2 norm of the distance between the two trajectories in phase space. This norm is computed using the distance at K + 1 discrete time points t 0 , , t K . Equation (5b) enforces the continuity of v across consecutive segments, Equation (5c) is the linearized version of Equation (1) around s, and Equation (5d) enforces the direction v to be normal to f . Finally, η is a time-dilation term between the reference and the shadowing trajectory. We briefly describe the solution of the minimization problem in Equation (5), and more details can be found in Reference [21,22]. We reformulate Equation (5) as
Minimize v i 1 2 i = 0 K v i 2 2 ,
subject to A v ̲ = b ̲ ,
where matrix A is
A = Φ 1 I Φ 2 I Φ K I v ̲ = v 0 v 1 v K b ̲ = b 1 b 2 b K ,
and Φ i + 1 = P t i + 1 ϕ t i , t i + 1 , b i + 1 = P t i + 1 t i t i + 1 ( ϕ τ , t ) f s d τ . Equation (8) is obtained by writing the state vector v ( t ) in t i t < t i + 1 as
v ( t ) = P t ϕ t , τ v ( t i ) + t i t ϕ t , τ f s d τ ,
where ϕ denotes the state transition matrix satisfying
d ϕ t , τ d τ = f u | τ ϕ t , τ ,
and P t is the projection operator which eliminates η from Equation (5)
P t = I f ( t ) f ( t ) T f ( t ) T f ( t ) .
Problem (7) is a standard minimization problem in linear algebra. It can be put into a Karush Kuhn Tucker (KKT) form:
I A T A 0 v ̲ w ̲ = 0 b ̲ .
Introducing the Schur complement S = A A T , the above system takes the form:
S w ̲ = Φ 1 Φ 1 T + I Φ 2 T Φ 2 Φ 2 Φ 2 T + I Φ 3 T Φ K Φ K Φ K T + I w 1 w 2 w K = b ̲ .
System (13) is solved with a Krylov solver, such as GMRES, by supplying matrix vector products S w ̲ . The system is usually ill-conditioned, and the Krylov solver requires a preconditioner to accelerate its convergence. In this work, we use the preconditioner developed in Reference [21], which results in
γ I + Γ S w ̲ = Γ b .
The introduction of matrix Γ and the regularization parameter γ reduces the condition number of S by orders of magnitude and accelerates the convergence of the iterative Krylov solver. After solving System (14), the sensitivity of G ¯ ( s ) is given by
d G ¯ ( s ) d s = 1 T i = 0 K 1 t i t i + 1 G u | t , v ̲ d t + 1 T i = 0 K 1 f i + 1 , v ̲ ( t i + 1 ) f i + 1 2 2 ( G ¯ G i + 1 ) + J ¯ s .
This approach to sensitivity analysis in chaotic systems has been shown to produce sensitivities that are representative of the physical properties of the underlying system and are in good agreement with finite differences. All gradients in this paper will be evaluated using the aforementioned algorithm.

3. Uncertainty Quantification with the gPC

In this section, we focus on the UQ methodology for chaotic systems described by Equation (1) for the QoI (2). Uncertainty is introduced through s which is modeled by a vector of m random independent variables ξ = ( ξ 1 , , ξ m ) , where each ξ i follows a given probability density function (PDF) w i ( ξ i ) , defined in the domain E i . Under these conditions, the QoI is also uncertain, becomes a function of ξ , i.e., G = G ( ξ ) , and we are interested in evaluating the statistical moments of G .
The vector ξ follows a PDF given by the product W = i = 1 m w i ( ξ i ) defined in the domain E = i = 1 m E i . The set of stochastic components ξ i defines a set of stochastic orthonormal bases ψ ( i ) = { ψ 0 ( i ) , ψ 1 ( i ) , } with a tensor product in the form Ψ : = i = 1 m Ψ ( i ) = { Ψ 0 , Ψ 1 , } that is orthogonal to the PDF W, i.e.,
Ψ j , Ψ k = E Ψ j Ψ k W d ξ = δ j k Ψ j , Ψ j ,
where δ j k denotes the Kronecker’s symbol, and repeated superscripts do not imply Einstein summation. Using this approach, an uncertain quantity F = F ( ξ ) can be represented in spectral form as
F ( ξ ) = i = 0 F i Ψ i ( ξ ) ,
which is commonly referred to as the PCE. This expansion is usually truncated after a finite number of terms, given by
q + 1 = ( m + c ) ! m ! c ! ,
where c is a user defined quantity, known as chaos order. The computation of the statistics of F requires computing the spectral coefficients F i . These coefficients can be computed through Galerkin projection, as in
F i = E F Ψ i d ξ .
In practice, the Galerkin integral is approximated by a numerical quadrature technique,
E F Ψ i d ξ j = 1 α ω j J ( ξ j ) Ψ i ( ξ j ) ,
that requires the evaluation of the QoI at points ξ j . The first 4 moments of F are computed through the following expressions:
Mean value E [ F ] = F 0 Standard deviation σ 2 [ F ] = E [ F 2 ] E 2 [ F ] Skewness S [ F ] = E [ F 3 ] 3 E [ F ] σ 2 [ F ] E 3 [ F ] σ 3 [ F ] Kurtosis K [ F ] = E [ F 4 ] 4 E [ F ] E [ F 3 ] + 6 E 2 [ F ] σ 2 [ F ] + 3 E 4 [ F ] σ 4 [ F ]
This framework can be used for conducting UQ either intrusively or non-intrusively. In this work, all Galerkin integrals are computed using Gauss Quadrature due to its numerical accuracy, although it is also common to use sparse grid approaches. For unsteady systems, that may or may not display chaotic behavior, the intrusive approach is not guaranteed to converge to correct statistics [24]. We illustrate this below with the aid of an example.
In the iPC method, the degrees of freedom are written as (refer to Equation (17)),
u ( ξ ) = i = 0 q u i Ψ i ( ξ )
and this expansion is inserted into the system equations as in
d d t ( i = 0 q u i Ψ i ( ξ ) ) = f ( i = 0 q u i Ψ i ( ξ ) , s ) .
To find the spectral coefficients of the expansion, a new set of equations, called the iPCE equations are derived by applying Galerkin projection to Equation (23),
d d t E ( i = 0 q u i Ψ i ( ξ ) ) Ψ i d ξ = E f ( i = 0 q u i Ψ i ( ξ ) , s ) Ψ i d ξ .
We now apply iPCE to the Lorenz attractor, a chaotic system commonly used as a test case,
d x d t = σ ( y x ) , d y d t = x ( ρ z ) y , d z d t = x y β z ,
where s = [ σ , ρ , β ] denotes the parameters vector (typical values are σ = 10 , ρ = 28 and β = 8 3 ). The system variables are written in their spectral form
[ x ( t ) , y ( t ) , z ( t ) ] = [ i = 0 q x i Ψ i ( ξ ) , i = 0 q y i Ψ i ( ξ ) , i = 0 q z i Ψ i ( ξ ) ] .
For demonstration purposes, we take m = 1 and c = 1 , and Equation (18) gives q = 1 . Uncertainty is introduced through β = β 0 + β 1 ξ . Substituting into the Lorenz system and applying Galerkin projections, we eventually get the system
d x 0 d t = σ ( y 0 x 0 ) d x 1 d t = σ ( y 1 x 1 ) d y 0 d t = x 0 ρ x 0 z 0 x 1 z 1 y 0 d y 1 d t = x 1 ρ x 0 z 1 x 1 z 0 y 1 d z 0 d t = x 0 y 0 + x 1 y 1 β 0 z 0 β 1 z 1 d z 1 d t = x 0 y 1 + x 1 y 0 β 0 z 1 β 1 z 0 .
This coupled system is integrated explicitly with a variable time-step Runge Kutta and unit initial conditions. We also performed MC simulations with N = 5000 samples over a time-frame T = 40 , which ensures the solution has become chaotic. The mean value (expectation) of the x ( t ) coordinate of the trajectory obtained from iPCE and Monte–Carlo (MC) are compared in Figure 1.
The MC result is obtained as follows. The Lorenz equations in (25) are integrated using N sample values of β that follow the given distribution. At each time instant t, the expectation of x ( t ) is computed from the N available samples and plotted as a function of time. Note that the expectation tends to 0 as t increases; this is due to the symmetry of the attractor. The iPC solution, however, starts to deviate substantially from the MC at around t = 15 . Increasing the chaos order c will not alleviate this deviation; it will simply delay it at later time.

4. Non-Intrusive PC in Chaotic Systems

In the niPC, instead of applying the PCE to the state variables, we apply it directly to the QoI, as in
G ¯ ( ξ ) = i = 0 G ¯ i Ψ i ( ξ ) ,
and the spectral coefficients are evaluated using a quadrature approach,
G ¯ i = E G ¯ ( ξ ) Ψ i d ξ = j = 1 α ω j G ¯ ( ξ j ) Ψ i ( ξ j ) .
The niPC requires evaluations of the QoI at the required quadrature points ξ j and does not suffer from the aforementioned problem of intrusive PC. This approach is therefore suitable to chaotic systems. It is easy to apply, and most importantly, time-averaged quantities G ¯ ( s ) are (usually) smooth functions of the parameters s and are therefore better suited for a spectral representation. This results in a smooth convergence of the spectral expansion and allows the evaluation of the statistics with low chaos orders c, which is important for keeping the computational cost low. More applications on niPC on chaotic systems will be discussed in later sections.
The same approach can be implemented to compute the statistics of the sensitivity of G ¯ ( s ) ,
d G ¯ d s i = 0 q d G ¯ d s i Ψ i ( ξ ) , d G ¯ d s i = j = 1 α ω j d G ¯ d s ( ξ j ) Ψ i ( ξ j ) .
The MSS algorithm, as described in Section 2, is employed in to evaluate d G ¯ d s ( ξ j ) at the quadrature points ξ j . In the following two sections, we apply niPC for time-averaged quantities and their sensitivities on two test cases and compare the results with MC simulations.

5. Application on the Lorenz-96 Model

We apply the aforementioned methodology in the Lorenz-96 equations [25], which model the evolution of atmospheric quantities, such as vorticity or temperature at a discrete periodic lattice representing a latitude circle on earth. The model is given by the system of N equations
d X i d t = ( X i + 1 X i 2 ) X i 1 X i + F ,
where i = 1 , , N . In this system, X i is a damping term, F is an external force, and the quadratic terms resemble the advection term and conserve kinetic energy. For more details on the solution of (31) and the treatment of X 1 , X 0 , and X N , the reader is referred to References [25,26]. This model is excellent for testing, since it is exhibits different dynamics depending on F; thus, it is ergodic, and its behavior displays similarities with turbulent fluid flows [26].
In Figure 2, we plot the maximum Lyapunov exponent of the system as a function of the forcing term F. Notice that the exponents are positive in the range of F considered, indicating that the Lorenz-96 model displays chaotic behavior. The exponents were computed using QR decomposition, as in Reference [27], and match well with the results from Reference [28].
The QoI considered is the time-averaged integral of the average value of X i , defined as
J ( F ) = 1 N T i = 1 N 0 T X i d t .
Uncertainty is introduced through F = 10 + 0.3 ξ , where ξ N ( 0,1 ) is Gaussian. To choose an appropriate chaos order c, we perform a spectral convergence analysis for the first four moments of Equation (32), which is shown in Figure 3. Notice that, for c = 4 , all four moments stabilize to the value predicted by a Monte–Carlo simulation with N = 10,000 samples.
In Figure 4, we compare the first four statistical moments of Equation (32) obtained using the Monte–Carlo and gPC for c = 4 . Notice that the Monte–Carlo is in good agreement with the gPC, indicating that this approach can be a valuable method for UQ in the time averages of chaotic systems. Note that the skewness and kurtosis are close to 0 and 3, respectively, as the values for Gaussian distribution. A possible explanation is that J ( F ) is the sum of random variables that follow the same distribution, and from the central limit theorem, it follows that the PDF of J ( F ) is Gaussian.
The same principle applies in the UQ of the sensitivity of QoI (32) with respect to F. This sensitivity is obtained using the MSS algorithm described in Section 2. To that end, Equation (31) is differentiated with respect to F
d V i d t = X i 1 V i + 1 + ( X i + 1 X i 2 ) V i 1 X i 1 V i 2 V i + 1 η ( X i + 1 X i 2 ) X i 1 X i + F ,
where V i = d X i d F . Equation (33) is used as a constraint in the MSS algorithm (refer to Equation (5c)). Results can seen in Figure 5, where the sensitivity obtained from MSS is compared with FD.
The convergence of the GMRes algorithm applied to system (14) is shown in Figure 6. The residuals are reduced by more than 8 orders of magnitude in 13 iterations, which demonstrates the effectiveness of the preconditioner.
By coupling the gPC with the MSS, the statistics of the sensitivity of the QoI can be computed efficiently. In Figure 7, the first four statistical moments of d J d F are computed and compared with the results from Monte–Carlo. The gPC is found to be in good agreement with the results from the Monte–Carlo simulation. Note that again the skewness and kurtosis are predicted to be close to the values of the Gaussian distribution.
The PDF computed by the gPC coupled with the MSS can be compared with the one produced by the Monte–Carlo where the sensitivities are evaluated using FD. The PDF predicted by the gPC is computed as in Reference [14]. This is seen in Figure 8, where the PDF from the Monte–Carlo has been approximated by connecting the midpoints of the bins of its histogram. Notice that the two PDFs are in good agreement.

6. Application on the van der Pol Oscillator

The UQ methodology is now applied to the van der Pol oscillator, given by
d 2 y d t 2 = y + β ( 1 y 2 ) d y d t .
This is a non–conservative oscillator, where the non-linearity is introduced through its damping coefficient. It was initially used to model limit cycle oscillations in electrical circuits, but it has been applied in other fields, as well, such us seismology and neuron modeling. The QoI is the L 8 norm of ω = d y d t , which behaves similarly to the L norm, and provides a metric for the magnitude of the peaks of ω . The objective function is written as
J 1 8 = lim T 1 T 0 T ω 8 d t 1 8 .
Equation (34) is converted into a set of two, coupled ODEs [29] and integrated from an initial condition at t = 50 until t = 0 , which ensures that the state variable u = [ y , ω ] has reached the attractor. Further integration then follows from 0 to T to obtain the reference trajectory. The attractor of the state u is seen in Figure 9.
Integrating Equation (34) forward in time is computationally inexpensive; therefore, we can easily compute the variation QoI with β for a wide range of values. This is shown in Figure 10; J 1 8 increases in magnitude as β increases.
The variation is smooth; therefore, it is expected that the gPC can be used to efficiently compute the statistics of J 1 8 . This is seen in Figure 11, where a comparison is made between the gPC and the MC. Note that the selection of the number of samples for the MC and the chaos order c was made through a spectral convergence study for the gPC and through a sample convergence study for the MC, similarly to the Lorenz-96 case. Uncertainty is introduced through β = β 0 + β 1 ξ , where β 0 varies in [ 0.2,2 ] with a step of δ b = 0.05 , ξ N ( 0 , 1 ) is Gaussian, and β 1 = 0.5 β 0 . Notice that the two methods are again in good agreement. It is worth mentioning that there is a slight deviation with respect to the value obtained in the absence of uncertainties (blue circles). This deviation, arising from a deterministic approach to chaotic systems, is important and results in overestimation of the QoI if uncertainties are not taken into account.
The smooth behavior shown in Figure 10 implies also a smooth variation of the sensitivity d J 1 8 / d β with respect to β . This is verified in Figure 12 (left), where the mean value is computed with FD for a time-length of T = 1000 . In the same figure, a comparison is made between the gPC coupled with MSS, and the Monte–Carlo for N = 5000 samples. Uncertainty conditions are the same as in Figure 11.
Notice that the gPC is in good agreement with the MC. A small deviation, similar to that found in Figure 11, is also observed in the sensitivities. Here, the deviation is found for even smaller uncertainties, which indicates an even stronger influence of the non-linearity present in the design space of Equation (35).
In Figure 13, we plot the reference trajectory and its shadowing counterpart, in phase space. The difference between the two trajectories is exaggerated by two orders of magnitude so that a visual can be produced. Notice that the two trajectories remain close, i.e., shadow, each other in phase space. It is exactly this property that allows us to compute reliable sensitivities in chaotic systems.

7. Discussion

The results from Section 5 and Section 6 indicate that the non-intrusive gPC is a valid method for quantifying uncertainties of time-averages and their sensitivities with respect to system parameters in chaotic systems. The gPC is an established method in propagating and evaluating uncertainties, yet the accuracy of the intrusive version deteriorates when applied to general unsteady (let alone chaotic) systems. This was demonstrated in the Lorenz attractor, where it was shown that the statistics of the trajectory cannot be captured correctly with an iPC approach. Approaches where the orthonormal basis is recalculated and the spectral order of the expansion is increased have appeared in the literature [30,31]; however, these approaches have high computational cost and increased complexity.
Non-intrusive approaches, on the other hand, do not suffer from the aforementioned problem. This is demonstrated on the Lorenz-96 system and the van der Pol oscillator, where the statistics of the time-averages are computed accurately for a low spectral chaos order c, the results being in good agreement with the MC. The convergence of the gPC for low c is due to the fact that the time-averages display smooth variation to system parameters, as well as due to ergodicity are independent of the initial conditions used. The computational cost of this method scales with the number of stochastic parameters. In such high-dimensional cases, correlations may exist between the uncertain parameters. These issues are dealt with in the gPC literature with the use of sparse grids, data–driven techniques, or variational approaches, [1,3,32].
Regarding the sensitivities, to properly evaluate their statistics, the gPC is coupled with the MSS algorithm, which can accurately compute the sensitivities of time-averages, even for a large number of uncertain parameters. The MSS algorithm is formulated in such a way to allow for the non-intrusive nature of the niPC to be preserved and for already developed software conducting the sensitivity analysis of the system to be employed. This coupling introduces a basic algorithmic framework for conducting UQ in chaotic systems and can exist as a basis for developing an approach to robust design in the presence of chaos.

Author Contributions

Conceptualization, K.D.K. and G.P.; software, K.D.K; validation, K.D.K and G.P.; writing—original draft preparation, K.D.K.; writing—review and editing, G.P.; supervision, G.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the President’s Scholarship Scheme at Imperial College London.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ahlfeld, R.; Ciampoli, F.; Pietropaoli, M.; Pepper, N.; Montomoli, F. Data-driven uncertainty quantification for Formula 1: Diffuser, wing tip and front wing variations. J. Autom. Eng. 2019, 233, 6. [Google Scholar] [CrossRef]
  2. Sasikumar, P.; Suresh, R.; Gupta, S. Stochastic model order reduction in uncertainty quantification of composite structures. Compos. Struct. 2015, 125, 21–34. [Google Scholar] [CrossRef]
  3. Pepper, N. Meta-modeling on detailed geography for accurate prediction of invasive alien species dispersal. Sci. Rep. 2019, 2019, 16237. [Google Scholar] [CrossRef] [PubMed]
  4. Knio, O.M.; Le Maître, O.P. Uncertainty propagation in CFD using polynomial chaos decomposition. Fluid Dyn. Res. 2006, 38, 616–640. [Google Scholar] [CrossRef]
  5. Chatzimanolakis, M.; Kantarakias, K.D.; Asouti, V.G.; Giannakoglou, K.C. A painless intrusive polynomial chaos method with RANS-based applications. Comput. Methods Appl. Mech. Eng. 2019, 348, 207–221. [Google Scholar] [CrossRef]
  6. Asmussen, S.; Glynn, P.W. Stochastic Simulation: Algorithms and Analysis. In Stochastic Modelling and Applied Probability; Springer: Cham, Switzerland, 2007. [Google Scholar]
  7. McKay, M. Latin hypercube sampling as a tool in uncertainty analysis of computer models. In Proceedings of the 24th Conference on Winter Simulation, Arlington, VA, USA, 13–16 December 1992. [Google Scholar]
  8. Morokoff, W.; Caflisch, R.E. Quasi-Monte Carlo integration. J. Comput. Phys. 1995, 122, 218–230. [Google Scholar] [CrossRef] [Green Version]
  9. Wiener, N. The homogeneous chaos. Am. J. Math. 1938, 60, 897–936. [Google Scholar] [CrossRef]
  10. Xiu, D.; Karniadakis, G.E. Modeling uncertainty in flow simulations via generalized polynomial chaos. J. Comput. Phys. 2003, 187, 137–167. [Google Scholar] [CrossRef]
  11. Le Maître, O.; Knio, O.M.; Najm, H.N.; Ghanem, R.G. A stochastic projection method for fluid flow I. basic formulation. J. Comput. Phys. 2001, 173, 481–511. [Google Scholar] [CrossRef] [Green Version]
  12. Salvetti, M.V.; Meldi, M.; Bruno, L.; Sagaut, P. Reliability of Large-Eddy Simulations: Benchmarking and Uncertainty Quantification. Direct Large-Eddy Simul. 2018, X24, 15–23. [Google Scholar]
  13. Blatman, G.; Sudret, B. Adaptive sparse polynomial chaos expansion based on least angle regression. J. Comput. Phys. 2011, 230, 2345–2367. [Google Scholar] [CrossRef]
  14. Oladyshkin, S.; Class, H.; Helmig, R.; Nowak, W. A concept for data driven uncertainty quantification and its application to a carbon dioxide storage in geological formations. Adv. Water Resour. 2011, 34, 1508–1518. [Google Scholar] [CrossRef]
  15. Lacor, C.; Dinescu, C.; Hirsch, C.; Smirnov, S. Implementation of Intrusive Polynomial Chaos in CFD Codes and Application to 3D Navier-Stokes; Springer: Cham, Switzerland, 2013. [Google Scholar]
  16. Fish, J.; Wu, W. A nonintrusive stochastic multiscale solver. Int. J. Numer. Methods Eng. 2011, 88, 862–879. [Google Scholar] [CrossRef]
  17. Luchini, P.; Bottaro, A. Adjoint Equations in Stability Analysis. Ann. Rev. Fluid Mech. 2013, 46, 493–517. [Google Scholar] [CrossRef]
  18. Lea, D.J.; Allen, M.R.; Haine, T.W.N. Sensitivity analysis of the climate of a chaotic system. Tellus A 2003, 52, 523–532. [Google Scholar] [CrossRef]
  19. Wang, Q.; Hu, R.; Blonigan, P. Least Squares Shadowing sensitivity analysis of chaotic limit cycle oscillations. J. Comput. Phys. 2011, 267, 210–224. [Google Scholar] [CrossRef] [Green Version]
  20. Palmer, K.J. Exponential Dichotomies, the Shadowing Lemma and Transversal Homoclinic Points. Dyn. Rep. 1988, 1, 265–306. [Google Scholar]
  21. Shawki, K.; Papadakis, G. A preconditioned Multiple Shooting Shadowing algorithm for the sensitivity analysis of chaotic systems. J. Comput. Phys. 2019, 398, 108861. [Google Scholar] [CrossRef] [Green Version]
  22. Blonigan, P.J.; Wang, Q. Multiple Shooting Shadowing for Sensitivity Analysis of Chaotic Systems and Turbulent fluid flows. In Proceedings of the 53rd AIAA Aerospace Sciences Meeting, Kissimmee, FL, USA, 5–9 January 2015. [Google Scholar] [CrossRef]
  23. Kantarakias, K.D.; Shawki, K.; Papadakis, G. Uncertainty quantification of sensitivities of time-average quantities in chaotic systems. Phys. Rev. E 2020, 101, 022223. [Google Scholar] [CrossRef] [Green Version]
  24. Branicki, M.; Majda, A.J. Fundamental limitations of polynomial chaos for uncertainty quantification in systems with intermittent instabilities. Commun. Math. Sci. 2013, 11, 55–103. [Google Scholar] [CrossRef]
  25. Lorenz, E.N.; Emanuel, K.A. Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model. J. Atmos. Sci. 1998, 55, 399–414. [Google Scholar] [CrossRef]
  26. Abramov, R.V.; Majda, A.J. Blended response algorithms for linear fluctuation-dissipation for complex nonlinear dynamical systems. Nonlinearity 2007, 20, 2793–2821. [Google Scholar] [CrossRef]
  27. Dieci, L.; Russell, R.D.; Van Vleck, E.S. On the Compuation of Lyapunov Exponents for Continuous Dynamical Systems. SIAM J. Numer. Anal. 1997, 34, 402–433. [Google Scholar] [CrossRef] [Green Version]
  28. Karimi, A.; Paul, M.R. Extensive chaos in the Lorenz-96 model. Chaos 2010, 20, 043105. [Google Scholar] [CrossRef] [Green Version]
  29. Boyd, S.P.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004. [Google Scholar]
  30. Choi, M.; Sapsis, T.P.; Karniadakis, G.E. A convergence study for SPDEs using combined Polynomial Chaos and Dynamically—Orthogonal Schemes. J. Comput. Phys. 2013, 245, 281–301. [Google Scholar] [CrossRef]
  31. Zygiridis, T.; Papadopoulos, A.; Kantartzis, N.; Antonopoulos, C.; Glytsis, E.N.; Tsiboukis, T.D. Intrusive polynomial-chaos approach for stochastic problems with axial symmetry. IET Microw. Antennas Propag. 2019, 13, 782–788. [Google Scholar] [CrossRef]
  32. Le Maître, O.P.; Knio, O.M. Spectral Methods for Uncertainty Quantification With Applications to Computational Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
Figure 1. Comparison of the mean value of x coordinate of the trajectory, as computed by intrusive Polynomial Chaos (iPC) and Monte–Carlo (MC) with N = 5000 samples. Uncertainty is introduced through β = β 0 + β 1 ξ , where ξ N ( 0 , 1 ) .
Figure 1. Comparison of the mean value of x coordinate of the trajectory, as computed by intrusive Polynomial Chaos (iPC) and Monte–Carlo (MC) with N = 5000 samples. Uncertainty is introduced through β = β 0 + β 1 ξ , where ξ N ( 0 , 1 ) .
Algorithms 13 00090 g001
Figure 2. Maximum Lyapunov exponent of the Lorenz-96 model for different values of the forcing term F (blue dots). The integration time is T = 500 and N = 40 . The red dots represent the exponents reported in Reference [28].
Figure 2. Maximum Lyapunov exponent of the Lorenz-96 model for different values of the forcing term F (blue dots). The integration time is T = 500 and N = 40 . The red dots represent the exponents reported in Reference [28].
Algorithms 13 00090 g002
Figure 3. Spectral convergence of the first four moments of the Quantity of Interest (QoI) (32) in the Lorenz-96 system. Here, T = 100 , F = 10 + 0.3 ξ , where ξ N ( 0 , 1 ) .
Figure 3. Spectral convergence of the first four moments of the Quantity of Interest (QoI) (32) in the Lorenz-96 system. Here, T = 100 , F = 10 + 0.3 ξ , where ξ N ( 0 , 1 ) .
Algorithms 13 00090 g003
Figure 4. Uncertainty Quantification (UQ) in the Lorenz-96 model ( N = 40 and T = 100 ), for the QoI (32). Comparison of statistical moments computed by Monte–Carlo and generalized Polynomial Chaos (gPC) for ( N = 40 and T = 100 ) c = 4 .
Figure 4. Uncertainty Quantification (UQ) in the Lorenz-96 model ( N = 40 and T = 100 ), for the QoI (32). Comparison of statistical moments computed by Monte–Carlo and generalized Polynomial Chaos (gPC) for ( N = 40 and T = 100 ) c = 4 .
Algorithms 13 00090 g004
Figure 5. Comparison of sensitivities as computed by MSS and FD. The MSS was conducted for a segment of size Δ T = 0.7 , K = 35 total segments and N = 40 .
Figure 5. Comparison of sensitivities as computed by MSS and FD. The MSS was conducted for a segment of size Δ T = 0.7 , K = 35 total segments and N = 40 .
Algorithms 13 00090 g005
Figure 6. Convergence of the Multiple-Shooting Shadowing (MSS) algorithm for the Lorenz-96 system. The MSS was conducted on K = 35 time segments on a trajectory of total length T = 50 , for N = 40 and F = 8 .
Figure 6. Convergence of the Multiple-Shooting Shadowing (MSS) algorithm for the Lorenz-96 system. The MSS was conducted on K = 35 time segments on a trajectory of total length T = 50 , for N = 40 and F = 8 .
Algorithms 13 00090 g006
Figure 7. UQ in the sensitivity of the QoI given by Equation (32) with respect to F. Comparison of statistical moments, as computed by Monte–Carlo coupled with Finite-Differences and gPC ( c = 4 ) coupled with the MSS algorithm. Uncertainty is introduced through F = 10 + 0.3 ξ where ξ N ( 0 , 1 ) . Here, N = 40 and T = 100 .
Figure 7. UQ in the sensitivity of the QoI given by Equation (32) with respect to F. Comparison of statistical moments, as computed by Monte–Carlo coupled with Finite-Differences and gPC ( c = 4 ) coupled with the MSS algorithm. Uncertainty is introduced through F = 10 + 0.3 ξ where ξ N ( 0 , 1 ) . Here, N = 40 and T = 100 .
Algorithms 13 00090 g007
Figure 8. Probability Density Function (PDF) comparison between Monte–Carlo ( N = 10000 samples) and gPC ( c = 4 ). Left: PDF of (32). Right: PDF of the sensitivity of Equation (32) w.r.t. F. The PDF from the Monte–Carlo is approximated by evaluating the histogram and connecting the bin midpoints. Here, T = 100 and N = 40 .
Figure 8. Probability Density Function (PDF) comparison between Monte–Carlo ( N = 10000 samples) and gPC ( c = 4 ). Left: PDF of (32). Right: PDF of the sensitivity of Equation (32) w.r.t. F. The PDF from the Monte–Carlo is approximated by evaluating the histogram and connecting the bin midpoints. Here, T = 100 and N = 40 .
Algorithms 13 00090 g008
Figure 9. Attractor of the state vector of Equation (34), for different values of β . Here, T = 100 .
Figure 9. Attractor of the state vector of Equation (34), for different values of β . Here, T = 100 .
Algorithms 13 00090 g009
Figure 10. Variation of the norm as in Equation (35) against β , for T = 500 .
Figure 10. Variation of the norm as in Equation (35) against β , for T = 500 .
Algorithms 13 00090 g010
Figure 11. Comparison of mean value (left) and standard deviation (right) of Equation (35), as computed by the gPC coupled with MSS and the Monte–Carlo for N = 7500 samples. In both cases, T = 100 . The chaos order is c = 4 .
Figure 11. Comparison of mean value (left) and standard deviation (right) of Equation (35), as computed by the gPC coupled with MSS and the Monte–Carlo for N = 7500 samples. In both cases, T = 100 . The chaos order is c = 4 .
Algorithms 13 00090 g011
Figure 12. Comparison of mean value (left) and standard deviation (right) of the sensitivity of (35), as computed by the gPC coupled with MSS and the Monte–Carlo for N = 7500 samples. In both cases, T = 100 . The chaos order is c = 4 .
Figure 12. Comparison of mean value (left) and standard deviation (right) of the sensitivity of (35), as computed by the gPC coupled with MSS and the Monte–Carlo for N = 7500 samples. In both cases, T = 100 . The chaos order is c = 4 .
Algorithms 13 00090 g012
Figure 13. Shadowing and reference trajectory. Here, β r e f = 0.3 , β s h a d = 0.301 , and T = 20 .
Figure 13. Shadowing and reference trajectory. Here, β r e f = 0.3 , β s h a d = 0.301 , and T = 20 .
Algorithms 13 00090 g013

Share and Cite

MDPI and ACS Style

Kantarakias, K.D.; Papadakis, G. Application of Generalized Polynomial Chaos for Quantification of Uncertainties of Time Averages and Their Sensitivities in Chaotic Systems. Algorithms 2020, 13, 90. https://doi.org/10.3390/a13040090

AMA Style

Kantarakias KD, Papadakis G. Application of Generalized Polynomial Chaos for Quantification of Uncertainties of Time Averages and Their Sensitivities in Chaotic Systems. Algorithms. 2020; 13(4):90. https://doi.org/10.3390/a13040090

Chicago/Turabian Style

Kantarakias, Kyriakos Dimitrios, and George Papadakis. 2020. "Application of Generalized Polynomial Chaos for Quantification of Uncertainties of Time Averages and Their Sensitivities in Chaotic Systems" Algorithms 13, no. 4: 90. https://doi.org/10.3390/a13040090

APA Style

Kantarakias, K. D., & Papadakis, G. (2020). Application of Generalized Polynomial Chaos for Quantification of Uncertainties of Time Averages and Their Sensitivities in Chaotic Systems. Algorithms, 13(4), 90. https://doi.org/10.3390/a13040090

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop