Stochastic Expansions Including Random Inputs on the Unit Circle

Stochastic Expansions Including Random Inputs on the Unit Circle

Brandon A. Jones111Assistant Professor, Department of Aerospace Engineering and Engineering Sciences, C0600, AIAA Senior Member. University of Texas at Austin, Austin, TX, 78712, USA    Marc Balducci222Graduate Research Assistant, Colorado Center for Astrodynamics Research, UCB 431. University of Colorado, Boulder, CO, 80309, USA

Stochastic expansion-based methods of uncertainty quantification, such as polynomial chaos and separated representations, require basis functions orthogonal with respect to the density of random inputs. Many modern engineering problems employ stochastic circular quantities, which are defined on the unit circle in the complex plane and characterized by probability density functions on this periodic domain. Hence, stochastic expansions with circular data require corresponding orthogonal polynomials on the unit circle to allow for their use in uncertainty quantification. Rogers-Szegő polynomials enable uncertainty quantification for random inputs described by the wrapped normal density. For the general case, this paper presents a framework for numerically generating orthogonal polynomials as a function of the distribution’s characteristic function and demonstrates their use with the von Mises density. The resulting stochastic expansions allow for estimating statistics describing the posterior density using the expansion coefficients. Results demonstrate the exponential convergence of these stochastic expansions and apply the proposed methods to propagating orbit-state uncertainty with equinoctial elements. The astrodynamics application of the theory improves robustness and accuracy when compared to approximating angular quantities as variables on the real line.






= equinoctial orbital elements
= stochastic ordinary differential operator
= semimajor axis (one of the equinoctial elements)
= space of complex numbers
= matrix comprised of expansion coefficients
= vector of polynomial chaos expansion coefficients
= number of random inputs
= expected value of
= imaginary part of complex number
= modified Bessel function of order
= imaginary number
= cost function
= measure of equatorial gravity perturbation
= Gaussian/normal distribution with mean and variance
= space of natural numbers
= number of surrogate training samples
= probability measure
= cardinality of
= maximum degree of terms in polynomial expansion
= parameter for Rogers-Szegő polynomials
= space of real numbers
= real part of complex number
= rank of separated representation
= weighting coefficients in separated representation
= space of complex numbers on unit circle
= Toeplitz matrix
= time
= vector of quantities of interest
= estimated value of variable
= matrix of propagated training samples
= von Mises distribution with location and concentration
= wrapped normal distribution with mean and variance
= multi-index in
= -dimensional hypercube - the image of random variables
= (discrete) Dirac delta function
= relative error/difference
= Verblunsky coefficient
= von Mises distribution concentration parameter
= mean longitude (one of the equinoctial elements)
= subset of multi-indices specifying tensor product of maximum degree and dimension
= mean of random variable
= gravitation parameter of the Earth
= vector of independent random variables
= joint probability measure
= Borel sigma algebra for probability space
= standard deviation of random variable
= characteristic function of order
= matrix of evaluated multi-variate orthogonal polynomials
= univariate polynomial of degree
= multi-variate polynomial defined by multi-index
= event space

I Introduction

Modern engineering problems require methods of tractable uncertainty quantification for problems with data, both inputs and outputs, in a variety of domains. For example, position is represented as a quantity on the real line. Attitude states may be parameterized via one of several methods (quaternions, Euler angles, etc.) that are defined over a domain that is periodic over one (or more) independent variables. The simplest form of these are circular data defined on the unit circle, and, hence, are periodic over (or an equivalent domain). Using methods of stochastic expansions for uncertainty quantification for problems that include circular data (or the higher-dimension analogues on the sphere, etc.) require an extension to existing theory to account for their differences when compared to data on the real line. This paper describes an approach to uncertainty quantification via stochastic expansions, namely Polynomial Chaos Expansions (PCEs) and Separated Representations (SR), when including inputs and outputs defined on the unit circle.

PCEs and SR provide a means for approximating the solution of a stochastic ordinary differential equation that is square-measurable, possibly non-Gaussian, with respect to the input uncertainties. Such techniques use a projection of the stochastic solution onto a basis of orthogonal polynomials in stochastic variables that is dense in the space of finite-variance random variables. To maintain efficiency in the expansion, the basis functions must be orthogonal with respect to the density function of the inputs. In cases of random variables that fit in the Wiener-Askey scheme, e.g., Hermite polynomials for Gaussian random inputs, selection of the basis is straightforward. However, previous research fails to identify the correct polynomial basis for random inputs on the unit circle.

A principle motivation of this work is in the context of astrodynamics problems. Orbit state uncertainty propagation is an active area of research, and the potential issues were first raised in junkins_1996 (). Proposed methods for representing and propagating uncertainty include state transition tensors (fujimoto_2012b, ), Gaussian mixtures (demars_2013a, ; horwood_2011, ), polynomial chaos (jones_2013a, ), separated representations (balducci_2017, ), differential algebra (valli_2013, ), and others. Applications of stochastic expansion-based methods to astrodynamics have been demonstrated for both Earth-centered and interplanetary missions (jones_2013a, ; jones_2013b, ; jones_2015b, ; balducci_2017, ; feldhacker_2016d, ). They also enable robust and reliable design via optimization under uncertainty (feldhacker_2016d, ). To date, most demonstrations of stochastic expansions to astrodynamics problems focus on Cartesian coordinates, and some work is required to rigorously apply them to uncertainty propagation with orbit element sets.

In junkins_1996 (), the authors demonstrated the use of classical orbital elements for improved uncertainty propagation, and the equinoctial elements have grown in popularity since they reduce complications with singularities. The various forms of the equinoctial elements include an angle describing location of the spacecraft in its orbit (e.g., see (hintz_2008, )). To account for circular data in the formulation of the state Probability Density Function (PDF) with the equinoctial elements, horwood_2014 () proposes a Gauss-von Mises density. To enable an equinoctial elements-based approach with stochastic expansions and, more generally, any problem with circular data, this paper describes how to formulate basis functions on the unit circle for the polynomial methods. Simulated results demonstrate the rapid convergence when using the appropriate orthogonal basis, and leverages the approach for orbit-state uncertainty propagation with equinoctial elements.

This paper presents a framework for including circular data and random inputs in stochastic expansions, and demonstrates their use for orbit-state uncertainty propagation. For circular data with probability density described by the Wrapped Normal Distribution (WND), the Roger-Szegő polynomial provide an orthogonal basis for use in stochastic expansions. Use of these polynomials requires recasting the random inputs into complex variables, and the paper discusses the implications of this in the context of surrogate approximations. For random inputs not described by the WND, a general framework is presented and numeric results use the proposed procedure with the circular random inputs characterized by the von Mises density. The paper includes a demonstration of rapid convergence (as a function of polynomial degree) of the mean and variance of the posterior PDF for a simple system. Results when propagating orbit-state uncertainty via equinoctial elements is then presented as a more practical application of the theory.

The paper outline is as follows. Section II states the problem formulation and describes the (general) use of stochastic expansions as a solution to the problem. Section III describes statistical descriptions of circular data and the orthogonal polynomials on the unit circle for use in stochastic expansions. Section IV describes their use in propagation of uncertainty via stochastic expansions, and Section V then demonstrates the performance of the approach. Finally, conclusions are discussed.

Ii Stochastic Expansions

ii.1 Problem Setup

Let be a probability space with event space , sigma algebra , and probability measure . Let defined on the probability space be random inputs with stochastic dimension . Assume that the elements of are independent but not necessarily identically distributed. Methods based on stochastic expansions seek to produce an approximate solution to the stochastic Ordinary Differential Equation (ODE) problem


where is the ODE operator, are the quantities of interest (e.g., propagated position and velocity), and is time. Random inputs may correspond to stochastic initial conditions and/or force model parameters.

This work considers the case where at least one element of is defined on a circle. Such quantities may be parameterize as angles, unit vectors, or by the circular variable where


i.e., numbers on the unit circle in the complex plane. Euler’s formula


where , relates the circular variable with the angle . Such circular data have a PDF defined on the unit circle and account for periodicity on the domain, e.g., . To leverage theory describing directional quantities on , the codomain of the expansions described in this paper are defined in the space of complex numbers. The following sections present a method for including random inputs corresponding to directional quantities in stochastic expansions. Note that the presentation of stochastic expansion methods is generalized to admit elements of defined as complex numbers to accommodate such Quantities Of Interest (QOIs). Note that we use as a generic angle variable for the sake of consistency between presented theory and future application to orbit-state propagation.

ii.2 Polynomial Chaos Expansions

A Polynomial Chaos Expansion (PCE) provides a means for generating an approximate solution to the stochastic ODE problem in (1) by projecting the solution onto a dense, multi-variate polynomial basis on the space of random inputs and orthogonal with regards to the density . This method was first proposed by wiener_1938 () for Gaussian random inputs, extended on by ghanem_2002 (); ghanem_1998 (); ghanem_1999a (); ghanem_1999 (), and later applied to a larger class of random inputs based on the Wiener-Askey scheme (xiu_2002a, ). Recent work focuses on extending such methods to a broader scope of problems with high dimension (doostan_2007b, ; doostan_2009, ; ma_2009, ; nouy_2010, ; doostan_2011, ; yang_2012, ; doostan_2013, ), computationally expensive simulation software (see discussion and motivation in narayan_2012 (); ng_2014 (); zhu_2014 (); narayan_2014 ()), nonlinear dynamics over long time spans (gerritsma_2010, ; wan_2005, ; wan_2006, ; jones_2013a, ; jones_2013b, ; jones_2015b, ), and/or discontinuities in the solution as a function of the random inputs (wan_2005, ; wan_2006, ; peng_2015, ).

A PCE approximates the solution for a quantity of interest vector via




is the maximum degree of the polynomial expansion, and the multi-variate basis functions are defined by the density (see Section II.4). For the sake of compact notation, dependence on time is hereafter removed. Elements of the set are multi-indices that define the multi-variate polynomial via


i.e., the multi-variate basis function is a tensor product of univariate basis functions each of degree and a function of a single random variable , and


is the number of terms in the PCE.

This work primarily leverages sampling-based (non-intrusive) methods for approximating the coefficients , which treats an existing, deterministic ODE solver as a black box. Let be the set of random input vectors, and the quantity of interest generated using evaluations of the prescribed ODE solver. The least-squares approach to solving for a PCE minimizes the least-squares cost function (e.g., see hosder_2006 ())


where denotes the conjugate transpose of . This yields the solution




Note that the matrix is only a function of the known random input set , and is the matrix of values produced via the deterministic ODE solver. For a least-squares solution, . While not leveraged in this work, methods based on compressive sampling allowing for generating a PCE with with the expansion is sparse doostan_2011 ().

ii.3 Separated Representations

The method of separated representations (SR), also known as canonical decompositions (CANDECOMP) or parallel factor analysis (PARAFAC), is similar to PCEs in that a solution to a stochastic ODE problem is approximated by projecting the solution onto a multi-variate polynomial basis. The formulation of SR, however, is different from a PCE. Instead of a sum along the multi-index, a surrogate based on SR decomposes a multi-variate function into a sum of products of univariate functions


where is the rank of the surrogate and are the univariate functions. The weighting coefficients are such that each has unit norm. These univariate functions, or factors, are composed via a sum of coefficients and, in the case of this paper, predetermined polynomials of an orthogonal basis. That is


where are unknown coefficients and form a polynomial basis as in the PCE formulation, but where denotes the degree of the univariate polynomial.

SR utilizes least squares regression in order to calculate the unknown coefficients. Unlike PCEs, the methodology of SR solves for one direction at a time, which reduces the overall computation algorithm to a series of linear optimization problems. In this alternating least squares (ALS) approach,


where the variable denotes the direction being solved for and . For the solution to (16), we seek to solve the normal equation


In the ALS process, for each direction, we solve for using (17), where the coefficients are organized as


the vector is the same as when generating a PCE but for only one QOI, and the matrix is represented in the block format


where is given by


and is the th component for sample . For a more thorough explanation of the ALS process, including a formulation for a vector, see beylkin_2009b () and balducci_2017 ().

ii.4 Selection of Basis Functions

Selection of the basis function influences convergence for a stochastic expansion. The solution that produces the fastest convergence, as a function of the number of terms in the expansion, leverages basis functions orthogonal with respect to the posterior distribution (e.g., see discussion in lemaitre_2010 () [pp. 35-36]). In general, this is not practical since the posterior PDF is not known a priori. Instead, basis functions orthogonal to the input probability measure are used, i.e.,


where denotes the complex conjugate of , and for the case of orthonormal polynomials. Typically, (21) is expressed assuming real numbers, but we generalize the property for use later in this paper. When this property is satisfied for the basis functions in (4), then


i.e., the approximation converges to in the mean-squares sense (cameron_1947, ). Basis functions for many common PDFs are well known as part of the Wiener-Askey scheme (xiu_2002a, ), e.g., Hermite and Legendre polynomials for Gaussian and uniform distributed random variables, respectively. To date, polynomials for PDFs common in the area of directional statistics have not yet been identified and applied to uncertainty quantification using stochastic expansions. When employing such orthogonal polynomials, both PCEs and SR admit analytic solutions for moments of the posterior distribution as a function of basis function inner products (e.g., see lemaitre_2010 (); doostan_2013 ()). The next sections discuss directional statistics and the appropriate set of orthogonal polynomials with regards to a given density.

Iii Random Inputs on the Unit Circle

iii.1 Directional Statistics

The field of directional statistics defines the fundamental theory to describe circular data. The subject began with studying distributions on the compass or in time, both of which may be modeled as functions on the unit circle. Over time, complexity increased to enable the study of probability and statistics in higher dimensions (e.g., see seminal papers fisher_1953 (); mardia_1975a (); bingham_1974 (); kent_1982 ()). This section presents pertinent properties of distributions on the unit circle and the probability densities of interest in the current work.

The random angle has a characteristic function


with the th trigonometric moments ()


where and denote the real and imaginary parts of , respectively. These are analogous to moments of the PDF seen for typical random variables and uniquely characterize any distribution on the circle (e.g., see (mardia_2000, , p. 26)). For empirically determined quantities from a collection of samples,


The mean direction and the circular standard deviation provide an analog to the related quantities on the real line and are functions of the first trigonometric moment. The circular mean of is


and the circular standard deviation is


Note, in this paper, there is no notational distinction between the circular and normal mean/standard deviation aside from the variable over which it is defined. The circular mean and standard deviation may be approximated empirically using Monte Carlo samples to estimate .

iii.2 Distributions on the Unit Circle

There are two analogs of the Gaussian distribution on the unit circle: the wrapped Normal density (mardia_2000, , pp. 50-51) and the von Mises distribution (vonmises_1918, ). The (one-dimensional) Wrapped Normal Density (WND)


is parameterized by the mean direction and the variance-like quantity . While similar to the familiar normal distribution on the real line, the sum over “wraps" the normal distribution around the unit circle and accounts for the angle rollover. The characteristic function of the wrapped normal distribution is


Initially proposed for applications in atomic physics, the von Mises distribution (VMD) is


where is also the mean direction, describes the concentration of the PDF, and is the modified Bessel function of the first kind and order . This PDF accounts for angle ambiguity through the cosine in the argument of the exponent. For the von Mises distribution

Fig. 1: Example von Mises and wrapped normal densities

Figure 1 presents the WND and VMD with and different values of and . The indicated value of is based on a nonlinear least squares fit to minimize differences in , . In the top case (small and large ), the von Mises case has slightly more density in the tails, whereas the wrapped normal is slightly wider in the primary lobe. The more concentrated case demonstrates the similarity between the two densities for small variance. The best use of each distribution varies, and pewsey_2005 () discusses the difficulty in discriminating between the two when given samples from one population. See collett_1981 () for a survey of comparisons between the two densities.

iii.3 Orthogonal Basis Functions

The following sections describe the polynomials to be used for stochastic expansions with random inputs in . The first section describes known polynomials that are orthogonal with respect to the WND and the following section describes generation of polynomials for any PDF with known characteristic moments.

iii.3.1 Orthogonal Polynomials on the Unit Circle

Orthogonal polynomials on the Unit Circle (OPUC) serve as the basis functions for stochastic expansions with random inputs . While this section provides a brief overview of OPUCs pertinent to this work, we refer the reader to simon_2004 () for a detailed introduction. The following sections describe specific OPUCs for a given density .

The Szegő recursion (e.g., see (simon_2004, , p. 56)) enables efficient and stable software implementation of OPUCs. Let denote an unnormalized form of the polynomial , and is a (unique) polynomial of degree that is orthogonal to with . Using the Szegő recursion,




and are dubbed the Verblunsky coefficients. These coefficients are unique for a given OPUC defined by measure in Eq. (21), and


For the form of the characteristic function in Eq. (23), the Toeplitz matrix for a given measure on is


As described in (szego_1975, , pp. 287-288),


which, when combined with Eq. (38) and computed for , allows for computing the Verblunsky coefficients. While this provides a mathematical framework for generating the OPUC for a given using the associated characteristic function, it can be sensitive to finite-precision arithmetic. This is especially true for measures highly concentrated on (see results and discussion in Section V), and are consistent with numeric generation of orthogonal polynomials on the real line (e.g., see gautschi_1982 ()).

iii.3.2 Rogers-Szegő Polynomials

Originally presented by Gabor Szegő (szego_1926, ) and based on the -Hermite polynomials of Leonard Rogers (rogers_1893, ; rogers_1894, ), the Rogers-Szegő polynomials are orthogonal with respect to . The (normalized) Rogers-Szegő polynomial of degree is (atakishiyev_1994b, )


with parameter and the -binomial coefficients


Note that different forms of these polynomials exist in the literature, e.g., see discussion in (simon_2004, , p. 87). The Verblunsky coefficients for the Rogers-Szegő polynomials are


These polynomials form an orthogonal basis with regards to the WND measure (e.g., see discussion in simon_2004 (); atakishiyev_1994b ()) given input argument


and parameter


Note that, via Euler’s formula, the Rogers-Szegő polynomials are trigonometric polynomials of the random angle . Figure 2 illustrates the real and imaginary components of the Rogers-Szegő polynomials for , , and . Note that the polynomials are periodic in .

Fig. 2: Rogers-Szegő polynomials of various degree and

Iv Stochastic Expansions with Random Inputs on the Unit Circle

For the case of , we leverage OPUC, which may be generated using a known . Note that, while this work only considers known analytic functions , approximate polynomials may be generated using approximate values computed via random samples from an unknown population. The following sections describe the implication of including angle quantities as random inputs or quantities of interest.

iv.1 Random Inputs and QOIs

This section outlines considerations when employing random inputs and/or QOIs on . To aid in understanding future results and clarifying sample generation, methods for simulating data on the unit circle are also included. Note that, for random angles, this presentation employs . With appropriate scaling, one may instead use to match other common random inputs. For ease in presentation, the former is employed in this paper.

Numeric results require simulating data from PDFs on the unit circle. Simulating samples with distribution requires generating samples from and mapping them to the unit circle. Let be a random sample of the standard Gaussian density . Then


This is indicative of the WND being a “wrapped" Gaussian density around the unit circle. Random inputs with the WND may be parameterized by or . Samples of the latter are related to the former via (46). When describing random samples via , then Hermite polynomials are the appropriate choice. For random inputs , then one should use the Roger-Szegő polynomials with input . Section V.2 presents similarities in the two approaches. Generation of samples based on the VMD does not allow for a simple transformation of standard random variables to quantities on the unit circle. Hence, random inputs must be directly generated for , and on the approach may be employed. A simple algorithm to generate samples may be found in best_1979 ().

The QOI may be parameterized by or the circular variable . Unless stated otherwise, all angular QOIs in this work are parameterized via . In such cases, or when using an OPUC for , then the coefficients for a given surrogate ( for the PCE and for SR) are in . The presentation of methods for generating the surrogate in previous sections is agnostic to real or complex QOIs, and accounting for the differences in software implementation is straightforward.

iv.2 Analytic Moments for PCEs

Upon generating a surrogate for system response as a function of the random inputs, select characteristics of the posterior density (e.g., mean and variance) are an analytic function of the expansion coefficients. For the case of QOI , the trigonometric moments, and thus the circular mean and standard deviation, are approximated via coefficients of the expansion. For the sake of simplicity, this presentation assumes a single QOI (i.e., not a vector ) and only random inputs on the unit circle.

Trigonometric moments for the case where QOI may be generated using the coefficients of a PCE. In the case of ,


since and for . To solve for the term of the characteristic function,


and the final equality results from orthonormal polynomials . Note that the equation for the second trigonometric polynomial differs from the second moment about the mean for QOIs on the real line, i.e., the sum is computed over all elements. Higher moment may be generated in a similar fashion, but will require the computation of polynomial triple products, etc.

Given the first trigonometric moment, then the PCE-estimated circular mean and STD are found via


While (55) and (56) assume a scalar QOI , these may be extended to vector quantities via element-wise operations.

iv.3 Analytic Moments for SR

As with PCEs, an SR surrogate is able to produce analytical expressions of the moments of the approximated function. To get the first moments, it is straightforward to show that


and the orthogonality of leaves the function solely dependent on the zeroth order coefficients and the normalizing constants. Similarly, when considering the second moment,


The trigonometric moments, circular mean, and circular standard deviation may be generated in a manner similar to the PCE given the SR solution for .

V Numerical Examples

The following numeric tests demonstrate the efficacy of using stochastic expansions with a random input on the unit circle. The first case demonstrates the exponential convergence expected when using the appropriate OPUC in a PCE. The remaining sections then demonstrate practical application of the approach by propagating orbit-state uncertainty with equinoctial elements.

v.1 Stochastic Expansion Convergence

This section demonstrates the exponential convergence achieved when using OPUC with random inputs defined on . Like the convergence rate demonstrations in xiu_2002a (), these tests consider the stochastic differential equation


with random decay-rate coefficient and deterministic solution


For this case, , where is a random variable on the unit circle, and time unit. The mean and variance of are


which, for a given density of , are computed numerically via quadrature using 1000 uniformly spaced nodes on the unit circle111See CIRCLE_RULE software found at [Accessed April 29, 2018]. This achieves accuracy in the baseline mean and variance approaching floating point error. Note that these are not the circular mean or standard deviation since is not an angle.

To assess performance of the OPUC basis while eliminating the selection of from influencing the result, PCE coefficients are numerically propagated forward to , and surrogate-derived statistics are compared to numeric approximations of Eqs. (65) and (66). Let and denote the PCE coefficients for and , respectively. In the following, multi-indices and follow the same mathematical definition as but use a different symbol to emphasize their use for difference PCEs. Direct propagation of the coefficients uses a Galerkin projection (see xiu_2002a () for details) to produce the system of linear differential equations




The initial conditions are


and the PCE for has coefficients


which is the analytic solution for a PCE such that . The ODE for the PCE coefficients (see Eq. (67)) may be used in any Runge-Kutta integrator with sufficient order and step size to achieve required accuracy. All cases use the common 4th order method with a step size of . The triple product values are pre-computed using quadrature integration with the appropriate polynomials and stored for use in propagating the PCE coefficients.

Fig. 3: Convergence of the PCE with random inputs from WND (top) and VMD (bottom).

Figure 3 illustrates the convergence of and as a function of expansion degree for four different distributions of . The WND and VMD for this case use the parameters identified in Figure 1 with “concentrated" and “diffuse" referring to larger and smaller values of , respectively. Error is quantified via


where denotes the quantity to characterize, i.e., mean or standard deviation. Baseline values of the moments for all cases (with enough digits of precision to illustrate differences) are provided in Table 1. Even the similar VMD and WND for the concentrated case yield non-trivial differences in moments. Results imply an exponential convergence in the expansion accuracy as increases, thereby demonstrating the efficacy of the OPUC for random angles.

Case Density
Diffuse WND
Concentrated WND
Table 1: Posterior Moments for intrusive PCE test cases
Fig. 4: Numeric stability of orthonormal polynomials for von Mises density.

Using the numeric procedure in Section III.3.1 to generate the OPUC can suffer from numeric issues for highly concentrated PDFs. Figure 4 demonstrates errors in the OPUC for the VMD quantified as the inner product as a function of and degree . Each evaluation of the integral (Eq. (21)) uses quadrature points. These inner products should equal zero, however finite precision arithmetic yields non-zero values for large . Numeric issues result from the condition number for the matrices in Eq. (40) when all . Developing a more numerically stable solution is designated as future work, but previous results for and still demonstrate solution convergence for that case. Note that polynomials with analytic solutions for the Verblunsky coefficients, such as the Rogers-Szegő polynomials, do not suffer from this issue.

v.2 Earth Orbit Cases

Of particular interest in this work is propagation of orbit-state uncertainty when using equinoctial elements. These elements include five quantities (, , , , ) defined on with the sixth being the mean longitude . Hence, propagation of equinoctial elements requires the presented framework with a mixture of QOIs in and . A series of test cases are presented, some of which focus on specific considerations when mixing polynomials and QOIs in different spaces.

Except when employing two-body only dynamics in Section V.2.1, all samples are propagated in time via special perturbations and TurboProp hill_2007 (). In such cases, equinoctial elements are converted to Cartesian position and velocity, propagated forward in time, and converted back to equinoctial elements to serve as the QOIs for a given surrogate. Propagation employs a Dormand-Prince 8(7) prince_1981 () embedded Runge-Kutta propagator with a relative tolerance of . The gravitational force model uses km/s and . For the sake of simplicity, transformation between the Earth-fixed and inertial coordinate systems is assumed to be a simple rotation about an inertially fixed rotation axis. While these scenarios are restricted to the main problem in Earth-centered orbit propagation for the sake of simplicity and duplication, previous demonstrations of surrogate methods for uncertainty propagation consider higher-fidelity dynamics jones_2013a (); jones_2015b (); balducci_2017 ().

Mean/Loc. 7444.0 km -7.07110 7.07110 7.07110 7.07110 33.59 deg
Variance 20 10 10 10 10 varies
Table 2: A priori PDF parameters for orbit test cases at

Table 2 summarizes the a priori mean and standard deviations employed in the following tests. These are loosely based on a test case in horwood_2011 (), but with eccentricity and inclination increased to and deg, respectively, to promote increased spatial variations in the gravitational field (due to ) over a single orbit. Uncertainty in the equinoctial elements remains unchanged when compared to horwood_2011 (), except where required for a given test. Initial PDFs for the first six elements are Gaussian with the provided standard deviation. With propagation times on the order of a day or more, this case yields a non-Gaussian posterior PDF for a initial deg. Hermite polynomials are used as the basis functions for all random inputs corresponding to the five elements defined on . The and any other deviations from Table 2 are described in the appropriate section. Note that any reference to a Root-Mean-Square (RMS) error in the following tests is based on a comparison of realizations of the QOI via evaluation of the propagator and the PCE when given the same random inputs.

v.2.1 Orbit State Propagation: -Only Case

To focus on propagation of uncertainty for , this section employs two-body only dynamics to propagate the angle:


with random inputs and corresponding to and , respectively. The initial state PDF is consistent with Table 2 with deg. Figure 5 depicts the propagated PDF for after 35.0 and 70.0 hours. Surrogates presented in this section use , and . For each surrogate method, three types of expansions are considered. Hermite polynomials and random inputs defined by are considered with QOIs and . The third combination considers Roger-Szegő polynomials with the circular variable . A comparable VMD would require . For reasons of numeric stability, OPUC for the VMD are not considered in this test case.

Fig. 5: Normalized histogram of propagated samples after 35 (top) and 70 (bottom) hours.
Method Basis QOI Rel. Error Mean Rel. Error STD RMS Error [deg]
PCE Hermite
SR Hermite
Table 3: Surrogate performance for -only case with hours

Table 3 presents the performance of the surrogates when approximating the PDF. The baseline circular mean and standard deviation using 10 Monte Carlo samples are approximately and , respectively. Relative error was defined previously in Eq. (71). For the surrogates leveraging the QOI , mean direction and the standard deviation may be generated using Eq. (55) and Eq. (56). Since no analytic solution (as a function of the expansion coefficients) is available for the mean direction and standard deviation given QOI , we use random samples of the surrogate (with the same random inputs used to generate the baseline solution) to compute similar values. Results demonstrate at least a 1000 improvement in the approximate mean and standard deviation when using , and greater improvement in the RMS error. The results for different basis functions (Hermite and Rogers-Szegő) with QOI exhibit only slight differences in the RMS error for the PCE surrogates. Otherwise, no significant difference is exhibited.

Fig. 6: PCE coefficients for the -only case with Rogers-Szegő and Hermite polynomials.

Figure 6 illustrates the effect that QOI selection has on the coefficients of the PCE. When using Rogers-Szegő polynomials with QOI , the PCE is converging as the number of terms increases. However, the Hermite polynomial with case is failing to converge, there by yielding the poor results seed in Table 3. Note that the Hermite with QOI is not included since it duplicates the Rogers-Szegő results. Using the Roger-Szegő polynomials demonstrates a clear dependence on a select few PCE terms. These correspond to terms with high degree for the random input corresponding to , and sensitivity to this input is expected (balducci_2017, ; horwood_2011, ). The unconverged PCE using the Hermite polynomial fails to provide a similar assessment of solution sensitivity.

Fig. 7: Surrogate accuracy over time for the -only case.

Figure 7 presents the performance of the surrogates as a function of time. Identifiers H and H refer to the basis (Hermite) and the QOI ( or ). The circular mean and standard deviation for the H surrogate uses 10 random samples of the PCE. Note that the difference in accuracy for surrogates with QOI is not visible to the scale of the figure (except as noted below), which is consistent with Table 3. Performance when using QOI reduces over time as the propagated PDF becomes more diffuse. Brief spikes in the H case depict those times where the PDF is split by the angle boundary (such as the case illustrated in the top image of Figure 5), and the duration of this reduction in error increases over time. At later times and for all surrogates, using a QOI instead of yields improved performance when compared to the H case. The spikes in relative error for when using QOI correspond to times where the mean direction is approximately zero, thereby making the relative error unstable. The RMS error of independent realizations of the surrogates further demonstrate the improved robustness and accuracy when using a more principled approach to surrogate design. While error increases over time, which is expected, it may be mitigated by increasing and . While a detailed comparison of PCE and SR is outside of the scope of this work, SR exhibits improved RMS accuracy for times less than approximately 15 hours.

v.2.2 Orbit State Propagation: SMA-Only Case

To demonstrate the effect of random inputs and polynomials in on QOIs in , this section focuses on the effects of uncertainty in () on . To allow for to influence over time, this case includes the perturbation in orbit propagation. The sensitivity of on will depend on inclination and eccentricity (see (vallado_2007, , pp. 653-654) for discussion), which motivated the change in mean state when compared to test cases in horwood_2011 (). Note that, for this case, the a priori PDF is the WND with deg, and the orbit is propagated for ten orbit periods (approx. 17.75 hours). All PCE solutions use samples with the Rogers-Szegő polynomials, and are compared to a Monte Carlo analysis with samples and a resulting precise to approximately four digits. Figure 8 presents the posterior (non-Gaussian) marginal PDF for using these Monte Carlo samples.

Fig. 8: Normalized histogram of propagated samples given random .
Fig. 9: Relative error for PCE-determined mean (top) and standard deviation (bottom).

Figure 9 presents the accuracy of the PCE-determined mean and standard deviation for as a function of . The Monte Carlo baseline refers to the empirically determined values, while the PCE baseline refers to the and from the PCE. The mean and standard deviation converge to approximately ten and four digits of agreement, respectively, when compared to the independent samples. This implies agreement between the solutions to the accuracy of the Monte Carlo analysis. The PCE baseline implies continued convergence to the true value as increases. While not presented in the interest of brevity, tests with samples demonstrate a agreement in for the PCE when and further reductions in differences when compared to the case.

Fig. 10: RMS error for PCE-determined as a function of .

Figure 10 demonstrates the RMS error as a function of . Since the PCE-produced is a complex number, this error is quantified three ways as designated in the legend. For this case, the imaginary and real components of error are comparable in magnitude and decay rapidly with . While not further explored in this work, using the imaginary component of error as a proxy for solution error is designated for future study. Normally, estimates of PCE accuracy require cross-validation with independent samples, and a method that does not require additional propagated samples is desired.

v.2.3 Full Orbit State Propagation

Propagation of uncertainty for the full equinoctial state is demonstrated in this section for two values of the a priori . All elements of the orbit state are considered stochastic, thus . The first case continues to use the relatively small uncertainty  deg with an a priori WND. To examine performance when the prior PDF for is a VMD, the second case employs a more diffuse density with . This avoids the numeric issues for large discussed with Figure 4. Using the same regression procedure used to identify similar densities in Figure 1, this scenario also employs an a priori WND with deg.

Fig. 11: Normalized histogram of propagated and components for the case with deg.

Figure 11 presents normalized histograms of and components of propagated samples after hours with the small initial angle uncertainty. For the sake of comparison, the figure includes the approximate posterior normal and wrapped normal densities for and , respectively, based on the empirical and . The posterior marginal density for is non-Gaussian, which results from short-period variations induced by the perturbation. Histograms are not provided for the other four equinoctial elements since the marginal PDFs are approximately Gaussian (to the scale of the resulting figure).

PCE coefficients with the Rogers-Szegő polynomials for (corresponding to random input ) are presented in Figure 12. For this case, and . A relatively small subset of PCE coefficients influence the propagated PDF, implying that may be decreased when combined with compressive sampling. This sparse PCE results from the statistical independence of the equinoctial elements, which was quantified in balducci_2017 () using SR. In all cases, the PCEs appear to have converged to at least five digits or more.

Fig. 12: PCE coefficients () for the case with deg.
Fig. 13: Normalized histogram of propagated and components for the diffuse case.
Fig. 14: PCE coefficients () for the diffuse case.

Figures 13 and 14 illustrate the posteriori marginal PDFs and PCE coefficients for the case with larger initial circular standard deviation. Given the more diffuse PDF and its sensitivity to the nonlinear dynamics, results for this case are presented after 24 hours. As seen in Figure 13, the propagated marginal PDF (approximated via samples) is non-Gaussian. PCEs are generated with and samples. Like the more concentrated case, the PCE coefficients imply a sparse approximation may be leveraged. Slight variations in the posterior solution are seen when comparing the PCE coefficients for the different prior densities for . While the degree of convergence varies with the QOI, all PCEs appear to be converged to 4 digits or more.

deg deg
Table 4: Relative error in for cases when compared to Monte Carlo

Table 4 quantifies the accuracy of the surrogate solutions for all cases. Assessment of the Monte Carlo solution indicates a convergence to approximately four digits of precision for with samples. As expected, the case with the smaller initial produces better agreement when compared to the Monte Carlo solution. Consistent with previous results, increases in and will demonstrate improved accuracy of the stochastic expansions for propagating uncertainty with equinoctial elements.

Vi Conclusions

Special considerations are required when using a surrogate method such as PCE or SR for uncertainty propagation with circular random inputs or quantities of interest. Using the circular variable to parameterize directional quantities of interest mitigates the latter issue, but polynomial-based methods require a basis orthogonal with respect to the input probability density functions. Random inputs described by the wrapped normal density employ the Rogers-Szegő polynomials, and a method for numerically generating polynomials allows for generalizing the approach to other densities that are sufficiently diffuse. Using these polynomials allows for computing directional statistics characterizing the posterior distribution as a function of the coefficients of the stochastic expansion. Expansions using these polynomials exhibit rapid convergence as a function of degree, and enable propagation of uncertainty in systems with angles as stochastic variables. This includes propagation of orbit-state uncertainty when using the equinoctial orbital elements.

Mr. Balducci’s work was funded by the NASA Science and Technology Research Fellowship program, contract NNX15AP41H. The authors also thank John Kent of the University of Leeds for his discussions on this work.



  • (1) Junkins, J. L., Akella, M. R., and Alfriend, K. T., “Non-Gaussian Error Propagation in Orbital Mechanics,” Journal of the Astronautical Sciences, Vol. 44, No. 4, 1996, pp. 541–563.
  • (2) Fujimoto, K., Scheeres, D. J., and Alfriend, K. T., “Analytical Nonlinear Propagation of Uncertainty in the Two-Body Problem,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 2, 2012, pp. 497–509.
  • (3) DeMars, K. J., Bishop, R. H., and Jah, M. K., ‘‘Entropy-Based Approach for Uncertainty Propagation of Nonlinear Dynamical Systems,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 4, 2013, pp. 1047–1057.
  • (4) Horwood, J. T., Aragon, N. D., and Poore, A. B., “Gaussian Sum Filters for Space Surveillance: Theory and Simulations,” Journal of Guidance, Control, and Dynamics, Vol. 34, No. 6, 2011, pp. 1839–1851.
  • (5) Jones, B. A., Doostan, A., and Born, G. H., “Nonlinear Propagation of Orbit Uncertainty Using Non-Intrusive Polynomial Chaos,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 2, 2013, pp. 430–444.
  • (6) Balducci, M., Jones, B. A., and Doostan, A., “Orbit uncertainty propagation and sensitivity analysis with separated representations,” Celestial Mechanics and Dynamical Astronomy, Vol. 129, No. 1-2, 2017, pp. 105–136.
  • (7) Valli, M., Armellin, R., Di Lizia, P., and Lavagna, M. R., “Nonlinear Mapping of Uncertainties in Celestial Mechanics,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 1, 2013, pp. 48–63.
  • (8) Jones, B. A. and Doostan, A., “Satellite Collision Probability Estimation Using Polynomial Chaos Expansions,” Advances in Space Research, Vol. 52, No. 11, 2013, pp. 1860–1875.
  • (9) Jones, B. A., Parrish, N., and Doostan, A., “Postmaneuver Collision Probability Estimation Using Sparse Polynomial Chaos Expansions,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 8, 2015, pp. 1425–1437.
  • (10) Feldhacker, J. D., Smith, J., Jones, B. A., and Doostan, A., “Multi-Element Trajectory Models for Satellite Tour Missions,” in “AIAA/AAS Astrodynamics Specialist Conference,” AIAA 2016-5502, Long Beach, California, 2016.
  • (11) Hintz, G. R., “Survey of Orbit Element Sets,” Journal of Guidance, Control, and Dynamics, Vol. 31, No. 3, 2008, pp. 785–790.
  • (12) Horwood, J. T. and Poore, A. B., “Gauss von Mises Distribution for Improved Uncertainty Realism in Space Situational Awareness,” SIAM/ASA Journal on Uncertainty Quantification, Vol. 2, No. 1, 2014, pp. 276–304.
  • (13) Wiener, N., “The Homogeneous Chaos,” American Journal of Mathematics, Vol. 60, No. 4, 1938, pp. 897–936.
  • (14) Ghanem, R. G. and Spanos, P. D., Stochastic Finite Elements: A Spectral Approach, Dover, New York, 2002.
  • (15) Ghanem, R. G. and Dham, S., “Stochastic Finite Element Analysis for Multiphase Flow in Heterogeneous Porous Media,” Transport in Porous Media, Vol. 32, No. 3, 1998, pp. 239–262.
  • (16) Ghanem, R. G. and Red-Horse, J., “Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach,” Physica D: Nonlinear Phenomena, Vol. 133, No. 1-4, 1999, pp. 137–144.
  • (17) Ghanem, R. G., “Ingredients for a General Purpose Stochastic Finite Elements Implementation,” Computer Methods in Applied Mechanics and Engineering, Vol. 168, No. 1-4, 1999, pp. 19–34.
  • (18) Xiu, D. and Karniadakis, G. E., “The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations,” SIAM Journal of Scientific Computing, Vol. 24, No. 2, 2002, pp. 619–644.
  • (19) Doostan, A., Iaccarino, G., and Etemadi, N., “A least-squares approximation of high-dimensional uncertain systems,” Tech. Rep. Annual Research Brief, Center for Turbulence Research, Stanford University, 2007.
  • (20) Doostan, A. and Iaccarino, G., “A least-squares approximation of partial differential equations with high-dimensional random inputs,” Journal of Computational Physics, Vol. 228, No. 12, 2009, pp. 4332–4345.
  • (21) Ma, X. and Zabaras, N., “An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations,” Journal of Computational Physics, Vol. 228, No. 8, 2009, pp. 3084–3113.
  • (22) Nouy, A., “Proper Generalized Decompositions and Separated Representations for the Numerical Solution of High Dimensional Stochastic Problems,” Archives of Computational Methods in Engineering, Vol. 17, No. 4, 2010, pp. 403–434.
  • (23) Doostan, A. and Owhadi, H., “A Non-adapted Sparse Approximation of PDEs with Stochastic Inputs,” Journal of Computational Physics, Vol. 230, No. 8, 2011, pp. 3015–3034.
  • (24) Yang, X., Choi, M., Lin, G., and Karniadakis, G. E., “Adaptive ANOVA decomposition of stochastic incompressible and compressible flows,” Journal of Computational Physics, Vol. 231, No. 4, 2012, pp. 1587–1614.
  • (25) Doostan, A., Validi, A., and Iaccarino, G., “Non-intrusive low-rank separated approximation of high-dimensional stochastic models,” Computation Methods in Applied Mechanical Engineering, Vol. 263, 2013, pp. 42–55.
  • (26) Narayan, A. and Xiu, D., “Stochastic Collocation Methods on Unstructured Grids in High Dimensions via Interpolation,” SIAM Journal of Scientific Computing, Vol. 34, No. 3, 2012, pp. A1729–A1752.
  • (27) Ng, L. W. T. and Willcox, K. E., “Multifidelity approaches for optimization under uncertainty,” International Journal for Numerical Methods in Engineering, Vol. 100, No. 10, 2014, pp. 746–772.
  • (28) Zhu, X., Narayan, A., and Xiu, D., “Computational Aspects of Stochastic Collocation with Multifidelity Models,” SIAM/ASA Journal on Uncertainty Quantification, Vol. 2, No. 1, 2014, pp. 444–463.
  • (29) Narayan, A., Gittelson, C., and Xiu, D., “A Stochastic Collocation Algorithm with Multifidelity Models,” SIAM Journal on Scientific Computing, Vol. 36, No. 2, 2014, pp. A495–A521.
  • (30) Gerritsma, M., van der Steen, J.-B., Vos, P., and Karniadakis, G., “Time-dependent generalized polynomial chaos,” Journal of Computational Physics, Vol. 229, No. 22, 2010, pp. 8333 – 8363.
  • (31) Wan, X. and Karniadakis, G. E., “An adaptive multi-element generalized polynomial chaos method for stochastic differential equations,” Journal of Computational Physics, Vol. 209, No. 2, 2005, pp. 617–642.
  • (32) Wan, X. and Karniadakis, G. E., “Multi-Element Generalized Polynomial Chaos for Arbitrary Probability Measures,” SIAM Journal on Scientific Computing, Vol. 28, No. 3, 2006, pp. 901–928.
  • (33) Peng, J., Uncertainty Quantification via Sparse Polynomial Chaos Expansion, Ph.D. thesis, University of Colorado Boulder, Boulder, CO, 2015.
  • (34) Hosder, S., Walters, R. W., and Perez, R., ‘‘A Non-Intrusive Polynomial Chaos Method for Uncertainty Propagation in CFD Simulations,” in “44th AIAA Aerospace Sciences Meeting and Exhibit,” AIAA 2006-891, Reno, Nevada, 2006.
  • (35) Beylkin, G., Garcke, J., and Mohlenkamp, M. J., “Multivariate Regression and Machine Learning with Sums of Separable Functions,” SIAM Journal of Scientific Computing, Vol. 31, No. 3, 2009, pp. 1840–1857.
  • (36) Le Maître, O. P. and Knio, O. M., Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics, Springer, 2010.
  • (37) Cameron, R. H. and Martin, W. T., “The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals,” Annals of Mathematics, Vol. 48, No. 2, 1947, pp. 385–392.
  • (38) Fisher, R., “Dispersion on a Sphere,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, Vol. 217, No. 1130, 1953, pp. 295–305.
  • (39) Mardia, K. V., “Statistics of Directional Data,” Journal of the Royal Statistical Society. Series B (Methodological), Vol. 37, No. 3, 1975, pp. 349–393.
  • (40) Bingham, C., “An Antipodally Symmetric Distribution on the Sphere,” The Annals of Statistics, Vol. 2, No. 6, 1974, pp. 1201–1225.
  • (41) Kent, J. T., “The Fisher-Bingham Distribution on the Sphere,” Journal of the Royal Statistical Society. Series B (Methodological), Vol. 44, No. 1, 1982, pp. 71–80.
  • (42) Mardia, K. V. and Jupp, P. E., Directional Statistics, John Wiley and Sons, Ltd., Chichester, England, 2000.
  • (43) von Mises, R., “Über die Ganzzahligkeit der Atomgewichte und Verwandte Fragen,” Physikalische Zeitschrift, Vol. 19, 1918, pp. 490–500.
  • (44) Pewsey, A. and Jones, M., “Discrimination between the von mises and wrapped normal distributions: Just how big does the sample size have to be?.” Statistics, Vol. 39, No. 2, 2005, pp. 81 – 89.
  • (45) Collett, D. and Lewis, T., “Discriminating Between the Von Mises and Wrapped Normal Distributions,” Australian Journal of Statistics, Vol. 23, No. 1, 1981, pp. 73–79.
  • (46) Simon, B., Orthogonal Polynomials on the Unit Circle, Part I: Classical Theory, American Mathematical Society, Providence, RI, 2004.
  • (47) Szegö, G., Orthogonal Polynomials, American Mathematical Society, Providence, RI, 4th ed., 1975.
  • (48) Gautschi, W., “On Generating Orthogonal Polynomials,” SIAM Journal on Scientific and Statistical Computing, Vol. 3, No. 3, 1982, pp. 289–317.
  • (49) Szegő, G., “Ein Beitrag zur Theorie der Thetafunktionen,” Sitzungsberichte der Preussischen Akademie der Wissenschaften, physikalisch-mathematische, pp. 242–252.
  • (50) Rogers, L. J., “Second Memoir on the Expansion of certain Infinite Products,” Proceedings of the London Mathematical Society, Vol. s1-25, No. 1, 1893, pp. 318–343.
  • (51) Rogers, L. J., “Third Memoir on the Expansion of certain Infinite Products,” Proceedings of the London Mathematical Society, Vol. s1-26, No. 1, 1894, pp. 15–32.
  • (52) Atakishiyev, N. M. and Nagiyev, S. M., “On the Rogers-Szegő polynomials,” Journal of Physics A: Mathematical and General, Vol. 27, No. 17, 1994, pp. L611–L615.
  • (53) Best, D. J. and Fisher, N. I., “Efficient Simulation of the von Mises Distribution,” Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 28, No. 2, 1979, pp. 152–157.
  • (54) Hill, K., TurboProp Version 3.2, Colorado Center for Astrodynamics Research, University of Colorado at Boulder, 2007.
  • (55) Prince, R. J. and Dormand, J. R., “High order embedded Runge-Kutta formulae,” Journal of Computational and Applied Mathematics, Vol. 7, No. 1, 1981, pp. 67–75.
  • (56) Vallado, D. A. and McClain, W. D., Fundamentals of Astrodynamics and Applications, Microcosm Press and Springer, Hawthorne, CA and New York, NY, 3rd ed., 2007.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description