Polynomial chaos based uncertainty quantification in Hamiltonian, multi-time scale, and chaotic systems

Uncertainty quantification techniques allow one to quantify output variability in the presence of parametric uncertainty. Typically, the moments of the output distributions are computed using sampling methods such as Monte Carlo [1], Quasi-Monte Carlo [2], and importance sampling [3]. Non-sampling approaches include response surface [4, 5] and polynomial chaos based methods [6]. Depending on the problem, different methods are applicable/appropriate in different scenarios. Polynomial chaos based techniques for propagating uncertainty have been used on a multitude of applications such as aeroelastic modeling [7], transport in heterogeneous media [8], Ising models [9], switching systems [10], combustion [11], fluid flow [12], and materials models [13], to name a few.

Here we study the properties and utility of using polynomial chaos expansions to propagate uncertainty through systems that have either Hamiltonian structure, multiple scales, or display chaos. We point out that polynomial chaos [6] and chaos theory [14] are unrelated areas. Originally proposed by Nobert Wiener [6] in 1938 (prior to the development of chaos theory—hence the unfortunate usage of the term chaos), polynomial chaos expansions are a popular method for propagating uncertainty through low dimensional systems with smooth dynamics. They rely on expanding random variables in terms of orthogonal basis functions [15]. Note that the orthogonal polynomials are chosen such that they are orthogonal to one another with respect to the prior distribution on the uncertain parameters [15]. For example, if the underlying distribution on the uncertain parameters is Gaussian, then the associated orthogonal polynomials are Hermite polynomials [16]. Similarly, if the underlying prior distribution is uniform, the associated orthogonal polynomials are Legendre [16]. In general, one can construct orthogonal polynomials for arbitrary distributions [15]. The advantage of polynomial chaos based techniques is that they provide exponential convergence for smooth processes with finite variance [17].

Chaos, on the other hand, refers to “aperiodic long-term behavior in deterministic systems that exhibits sensitive dependence on initial conditions” [14]. Chaos theory has been applied to a wide variety of applications such as fluid turbulence [18], celestial dynamics [19], and weather modeling [20]. It is important to point out that although the dynamics has sensitive dependence on initial conditions, it is inherently deterministic. In other words, no associated parametric uncertainty is required to observe chaos.

In this work, we present three new results. In the first part, we show that the dynamical systems that one gets by applying polynomial chaos expansions to Hamiltonian systems with uncertain parameters are also Hamiltonian. To do this, we first perform a polynomial chaos expansion of the generalized coordinates and conjugate momenta and find the evolution equations for the coefficients. We consider the expansion coefficients of the generalized coordinates as a new, larger set of “uncertain” generalized coordinates. By considering the averaged Hamiltonian (over parameter space) as a function of the expansion coefficients, we show that, for each of these new generalized coordinates, the coefficient in the expansion of the corresponding conjugate momentum and to the corresponding order is itself the conjugate momentum relative to the average Hamiltonian, thus demonstrating that the Hamiltonian structure in the derived differential equations is preserved. We then connect this result with the volume-preserving property of Hamiltonian flows to show that any finite polynomial chaos expansion of a harmonic oscillator with uncertain frequency must fail at long times, regardless of how many terms are kept. In the second part, we demonstrate the application of polynomial chaos to systems with multiple time scales using perturbation theory [21]. We demonstrate how the uncertain parameters influence the averaged dynamics of the dynamical system. In particular, we show that uncertainty can be propagated through the averaged equations instead of through the original equations, thus avoiding the computational burden of simulating a stiff system. In the third part, we apply polynomial chaos to a chaotic dynamical system (forced Duffing oscillator) [22] and demonstrate that the resulting equations for the coefficients are also chaotic. We then show that chaotic dynamics significantly reduce the efficacy of polynomial chaos expansions at propagating uncertainty.

Starting with a complete probability space $\Gamma$ given by $(\Omega,\mathcal{F},\mathbb{P})$, where $\Omega$ is the sample space, $\mathcal{F}$ is the $\sigma$-algebra on $\Omega$ and $\mathbb{P}$ is a probability measure, let $L_{2}(\Gamma,X)$ denote the Hilbert space of square-integrable, $\mathcal{F}$-measurable, $X$-valued random elements. Then one can, in general, define a polynomial chaos basis $\{\psi_{k}(\lambda(\omega))\}$, where $\lambda(\omega)$ is a random vector, $\omega\in\Omega$, and $k=(k_{1},k_{2},\dots)$ is a vector of non-negative indices. We denote the probability density function of the random vector $\lambda$ by $\rho(\lambda)$.

Generalized polynomial chaos (gPC) [15] provides a framework for representing second-order stochastic processes $\kappa\in L_{2}(\Gamma,X)$ for arbitrary distributions of $\lambda$ by the following expansion:

where $|k|=\sum_{i}k_{i}$ is the sum of the indices of $k$ and $\psi_{k}(\lambda)$ are orthonormal polynomials on $\Gamma$ with respect to $\rho(\lambda)$. Restricting our formalism to $\mathbb{R}^{n}$ (relevant for this work) we get the orthonormality is given by

where $\delta_{ik}$ is the Kronecker delta product. Depending on $\rho(\lambda)$ one can generate an appropriate orthogonal basis for representing $\kappa(\lambda)$. As mentioned earlier, if $\rho$ is Gaussian, then the appropriate polynomial chaos basis is the set of Hermite polynomials; if $\rho$ is the uniform distribution, then the basis is the set of Legendre polynomials. For details on the correspondence between distributions and polynomials see [23, 24]. A framework to generate polynomials for arbitrary distributions has been developed in [15]. The advantage of using polynomial chaos is that it provides exponential convergence in smooth processes [17]. However, the approach suffers from the curse of dimensionality, making them infeasible for problems with more than a handful of parameters. To mitigate the curse of dimensionality, sparse grid techniques [25, 26, 27], iterative methods [28, 29, 30], regression based algorithms [31, 32], hierarchical methods [33], dimensionality reduction based techniques [34, 35] have been developed.

In practice, the expansion in Eqn. 1 is truncated at a particular order, say, $r$. One can then use Galerkin projections to obtain a set of differential equations for the coefficients $a_{k}$ in Eqn. 1 [15]. Typically, low order truncations are found to capture the uncertainty in smooth systems [17].

Consider a system described by the Hamiltonian $H(q,p;\lambda)$, where $\lambda$ is a vector of uncertain parameters with probability density $\rho(\lambda)$. The generalized coordinates and momenta $q_{i}$ and $p_{i}$ ($i=1,\ldots,N$) satisfy Hamilton’s equations.

The generalized polynomial chaos (gPC) expansion of the coordinates and momenta is,

where the $\psi_{k}$ form an orthonormal basis with respect to the density $\rho$ (see Eq. 2).

The gPC coefficients $Q_{ik}$ and $P_{ik}$ follow deterministic equations, obtained by projecting the equations of motion along $\psi_{s}$

Inserting the gPC expansions and using the orthonormality condition (2) we obtain,

Let us define the average Hamiltonian $\hat{H}$

By using the gPC expansion of $q$ and $p$ we can consider $\hat{H}$ as a function of $Q$ and $P$, where $Q$ and $P$ denote the sets of coefficients $Q_{ik}$ and $P_{ik}$, respectively.

Note that the proof depends only on the form of the expansion and does not require that the expansion be complete. In other words, the coefficients of a truncated expansion will also form a (finite) Hamiltonian system relative to the average Hamiltonian when expressed as a function of the truncated expansion coefficients. Hence, polynomial chaos expansions when applied to Hamiltonian systems are also Hamiltonian. This result is not only interesting but also has practical implications. In particular, if the underlying system is Hamiltonian and one desires to propagate uncertainty using polynomial chaos, symplectic integrators [36] will be needed to maintain numerical accuracy for long times.

We now illustrate the preservation of Hamiltonian structure on the Duffing oscillator with parametric uncertainty.

Systems with multiple time scales are prevalent in a wide variety of applications related to smart grids [38, 39, 40], building systems [28], and micromechanical oscillators [41, 42, 43], to name a few. Simulating these systems is challenging due to their inherent stiffness [44]. The method of multiple scales (or averaging) is a very popular approach for simulating such systems. The approach typically involves perturbing off a dynamical system whose solution can be computed in closed form [21]. Note that this approach is applicable only in the scenario that the perturbation is small $O(\epsilon)$. The method of multiple scales captures the dynamics of the system on an $n-1$ dimensional section transversal to the flow, also known as the Poincaré section [22]. For a detailed discussion on the method of multiple scales or averaging theory, see [21, 22].

To the best of our knowledge, no attempt has been made to extend polynomial chaos based methods to systems with multiple time scales using the method of multiple scales. Here we apply polynomial chaos to the two-time scale system given below,

where $\epsilon\delta$ is the system damping, $1$ and $\epsilon\beta$ are the linear and nonlinear stiffnesses respectively, and $\epsilon\gamma$ and $\omega$ are the forcing amplitude and frequency respectively. Note that in the above system, we assume that $\epsilon$ is a small parameter (i.e. $\epsilon\ll 1$). We assume that $\gamma=\gamma_{0}+\sigma(\gamma)\eta$ is an uncertain parameter , where $\gamma_{0}=1.0$ is the mean of $\gamma$ and $\sigma(\gamma)=0.1$ is its standard deviation ($\eta$ is a normal random variable with zero mean and unit variance).

Using the two time scales as $\xi=\omega t$ and $\chi=\epsilon t$, one can derive the averaged equations for the system [21, 22]. This is done by substituting $\frac{d}{dt}=\omega\frac{\partial}{\partial\xi}+\epsilon\frac{\partial}{\partial\chi}$, $\frac{d^{2}}{dt^{2}}=\omega^{2}\frac{\partial^{2}}{\partial\xi^{2}}+2\omega\epsilon\frac{\partial^{2}}{\partial\xi\partial\chi}+\epsilon^{2}\frac{\partial^{2}}{\partial\chi^{2}}$, and $q(\xi,\chi)=q_{0}(\xi,\chi)+\epsilon q_{1}(\xi,\chi)+\ldots$ in Eqn. 24. Collecting terms, we obtain

The solution to Eqn. 25, is $q_{0}(\xi,\chi)=A(\chi)\cos\xi+B(\chi)\sin\xi$. Substituting the solution into Eqn. 26 and imposing that there be no secular terms [21, 22] yields the averaged equations

Note that the above dynamical system captures the dynamics on the Poincaré section of the original system [22]. From here on we take $\beta=1$. We will also focus on the deterministic initial condition $q=2$, $\dot{q}=0$, so $A(0)=2$ and $B(0)=0$.

We now apply a polynomial chaos expansion to Eq. 27, to first order:

where $H_{0}(\eta)=1$ and $H_{1}(\eta)=\eta$ are the first two probabilist’s Hermite polynomials. Substituting Eq. 28 into Eq. 27 gives

with initial condition $(a_{0},b_{0},a_{1},b_{1})=(2,0,0,0)$.

For comparison purposes, we also do an equivalent polynomial chaos expansion of the original two-time equations, defining $x=q$, $y=\dot{q}$, and doing a first order expansion

The resulting system is

with initial condition $(x_{0},y_{0},x_{1},y_{1})=(2,0,0,0)$. In what follows, we choose $\omega=1$ and, more importantly, $\delta=0$. This is done in order to avoid having a system in which the dissipation artificially reduces the overall error.

Figure 2 shows the error in mean and standard deviation of both PC expansions compared with Monte Carlo simulations with $10^{3}$ samples of the original two-time system, evaluated at the Poincaré sections where the forcing is maximal: $t=2\pi n$ ($n=0,1,2,\ldots$). Solutions of both PC expansions and the Monte Carlo trajectories of the original system were obtained using Matlab’s ode45 solver with relative tolerance of $10^{-6}$. The error has two sources: the time averaging and the truncation of the polynomial chaos expansion. As we decrease $\epsilon$, the error from time averaging decreases and the main source of error becomes the truncation of the polynomial expansion.

As $\epsilon$ decreases, the two expansions yield increasingly similar results, but the PC expansion of the original equation becomes more expensive to compute, scaling as $1/\epsilon$, because the solver needs to trace each fast oscillation, even if we’re only interested in the slow evolution of the Poincaré section. Figure 3 shows the number of function evaluations required by this expansion as $\epsilon$ decreases (solid line) compared with the $\epsilon$-independent behavior of the averaged PC.

We now demonstrate that dynamics of the coefficients of the polynomial chaos expansions can be chaotic, if the underlying dynamical system is chaotic. We will then show that if the underlying system is chaotic, the applicability of polynomial chaos is significantly reduced. While this is a natural result, it is not obvious. PC aims to capture the moments of the output distribution and not individual trajectories of the underlying dynamical system. Since mixing introduces averaging that could, in principle, smooth out the moments and allow them to be captured accurately, it is not obvious that chaos would necessarily make predictions worse.

For this demonstration, we pick the forced Duffing oscillator [22] given by

where $\delta=0.2$, $\gamma=0.3$, $\omega=1.0$, and $\lambda=-1.0$. Note that the above equation (Eq. 31) is the same as Eq. 5 with the addition of damping and forcing terms. We can write Eq. 31 in the form

The dynamics of the above system have been studied extensively (see, e.g., [22]). For the forced Duffing oscillator, the Poincaré section is given by taking “snapshots” of the system at phase $\phi=0$, where $\phi=(\omega t\,\mbox{mod}\,2\pi)$. The intersection of a single trajectory with the Poincaré section can be seen in Fig. 4, starting from the initial condition $(q,p)=(1,0)$.

The dynamics of the forced Duffing oscillator (at the parameter values given above) is well known to be chaotic [22]. In fact, one can numerically compute Lyapunov exponents ($\Xi$) [14] for the above system and show that they are positive. Note that $\Xi>0$ is considered to be the signature of a chaotic system since it implies that the system response is sensitive to initial conditions. We find that the nominal system gives $\Xi\approx 0.93$, hence (numerically) implying the existence of chaos.

Let us now assume that $\lambda$ is normally distributed. Let $\lambda=\lambda_{0}+\sigma\eta$, where $\lambda_{0}=-1.0$ is the mean of $\lambda$ and $\sigma=0.1$ is its standard deviation. Is is easy to see that $\eta$ will now become a normally distributed random variable with zero mean and unit standard deviation. Since $\eta$ is normally distributed, we use Hermite polynomials in our expansion [6]. In Eq. 32 we use the expansion in Eq. 8. Truncating the expansion at $r=1$ gives the following set of differential equations:

Note that there is nothing special about order $r=1$, and the same procedure can be repeated for any $r$. The initial condition on the generalized coordinates $q$ and conjugate momenta $p$ gets incorporated into the initial conditions on the coefficients of expansion: $Q_{i}$ and $P_{i}$. The Poincaré section for Eqs. 33 are shown in Fig. 5. In this case, the stretching and folding structure of the Duffing oscillator is not as evident as in Fig. 4. However, the resulting dynamical system in Eqs. 33 has a Lyapunov exponent of $\Xi\approx 0.73$, implying the persistence of sensitive dependence to initial conditions. Note that the route to chaos [22] for the the original Duffing oscillator is well known. In [45, 46], the authors numerically demonstrate that the forced Duffing oscillator becomes chaotic due to a sequence of period doubling bifurcations. Due to the onset of chaos, the solution becomes increasingly difficult for polynomial chaos to capture. We point out that polynomial chaos is known to suffer from an inability to track output distributions for long term simulations [47].

The reason for the inability of polynomial chaos to track the output distribution lies in the increasingly oscillatory nature of the solution $q(t;\lambda)$ in terms of the uncertain parameter $\lambda$. In other words, any finite expansion in Eqn. 8 will fail at some $t$, since $q(t,\lambda)$ is too oscillatory in terms of $\lambda$. The greater the oscillatory nature of the output in terms of $\lambda$, the worse polynomial chaos performs. In [10], the oscillatory nature of the output is again found to adversely impact the propagation of uncertainty through hybrid dynamical systems. However, we find that chaotic dynamics exacerbates this phenomenon. In particular, due to the coexistence of periodic orbits of different periods along with the chaotic attractor, the solution is found to rapidly become oscillatory with respect to $\lambda$ (depicted in Fig. 6). Polynomial chaos is unable to track the first moment (mean) of $q(t;\lambda)$ beyond certain time (10 secs in Fig. 7, 25 secs in Fig. 9). In contrast to the forced Duffing oscillator, polynomial chaos is able to accurately track the mean of $q$ in the undamped and unforced Duffing oscillator with an order of expansion of $r=1$ (see Figs. 8 and 10). Note that all parameters and initial conditions are held constant here (except for the removal of the forcing and damping terms). Hence, one needs to be careful when applying polynomial chaos to systems that are chaotic. We point to a caveat that if the initial condition is chosen close to the stable and unstable manifolds of the saddle equilibrium $(0,0)$, polynomial chaos performs poorly on the undamped, unforced oscillator case due to the discontinuity associated with the basin boundary [10].

When propagating uncertainty through chaotic systems with uncertain initial conditions, polynomial chaos again will need to be used carefully. Assume that $\lambda$ is not uncertain anymore, but instead the initial conditions are normally distributed as $(q,p)=(1,0)+(\sigma\eta,0)$, where $\eta$ is a Gaussian variable with zero mean and unit variance. The first order expansion yields the following system:

Note that the uncertainty in initial conditions does not appear explicitly in these equations, but rather enters through the initial condition $Q_{1}(0)=\sigma$. For the purpose of our simulations we take $\sigma=0.1$. The resulting Poincaré sections are shown in Fig. 11. The Lyapunov exponent is numerically found to be $\approx 0.85$, suggesting the persistence of chaos in the resulting polynomial chaos equations. This implies that any long term simulation that aims to track the output distribution will also suffer from problems of round-off in the initial conditions (given that the distribution on the initial condition will require computation of the initial conditions of the coefficients).

Polynomial chaos is slowly becoming an established and popular approach for propagating uncertainty through smooth systems. Every year researchers use the approach to propagate uncertainty through a wide variety of engineering [7, 8, 11, 13] and biological systems [48]. A systematic study on the properties and applicability of polynomial chaos to systems based on their structure and dynamics appears to be lacking.

In this work, we presented three main results. In the first part, we proved that when polynomial chaos is applied to Hamiltonian systems, the resulting equations are also Hamiltonian, even when the expansion is truncated. This is important, as it implies that structure in Hamiltonian systems is not only inherited by the new equations but also require the use of structure-preserving integrators [36] to accurately propagate uncertainty. We also used the volume-preserving property of Hamiltonian systems to show that, on a particular example, a finite expansion must fail at long times, regardless of the order of the expansion. In the second part, we show that polynomial chaos may be applied to the averaged equations of a forced two-time system, allowing much faster uncertainty propagation than polynomial chaos on the original system. As the time scale separation increases, both the computational advantage as well as the quality of the approximation improves. In the third part, we demonstrate that polynomial chaos also inherits chaotic dynamics from underlying systems. The presence of chaos is shown to negatively influence the applicability of polynomial chaos. It reduces the length of time that polynomial chaos accurately tracks the output distributions and complicates computations when there is uncertainty in initial conditions.

The authors would like to thank Amit Surana, as well as the anonymous reviewers for constructive comments.