Stochastic Expansions Including Random Inputs on the Unit Circle
Abstract
Stochastic expansionbased 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. RogersSzegő 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 orbitstate 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.
\@ifdim
>
Nomenclature
=  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 RogersSzegő 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  
=  multiindex 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 multiindices 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 multivariate orthogonal polynomials  
=  univariate polynomial of degree  
=  multivariate polynomial defined by multiindex  
=  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 higherdimension 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 squaremeasurable, possibly nonGaussian, 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 finitevariance 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 WienerAskey 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 expansionbased methods to astrodynamics have been demonstrated for both Earthcentered 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 Gaussvon Mises density. To enable an equinoctial elementsbased 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 orbitstate 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 orbitstate uncertainty propagation. For circular data with probability density described by the Wrapped Normal Distribution (WND), the RogerSzegő 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 orbitstate 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
(1) 
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
(2) 
i.e., numbers on the unit circle in the complex plane. Euler’s formula
(3) 
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 orbitstate 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, multivariate 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 WienerAskey 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
(4) 
where
(5)  
(6) 
is the maximum degree of the polynomial expansion, and the multivariate 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 multiindices that define the multivariate polynomial via
(7) 
i.e., the multivariate basis function is a tensor product of univariate basis functions each of degree and a function of a single random variable , and
(8) 
is the number of terms in the PCE.
This work primarily leverages samplingbased (nonintrusive) 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 leastsquares approach to solving for a PCE minimizes the leastsquares cost function (e.g., see hosder_2006 ())
(9) 
where denotes the conjugate transpose of . This yields the solution
(10) 
where
(11)  
(12)  
(13) 
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 leastsquares 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 multivariate polynomial basis. The formulation of SR, however, is different from a PCE. Instead of a sum along the multiindex, a surrogate based on SR decomposes a multivariate function into a sum of products of univariate functions
(14) 
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
(15) 
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,
(16) 
where the variable denotes the direction being solved for and . For the solution to (16), we seek to solve the normal equation
(17) 
In the ALS process, for each direction, we solve for using (17), where the coefficients are organized as
(18) 
the vector is the same as when generating a PCE but for only one QOI, and the matrix is represented in the block format
(19) 
where is given by
(20) 
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. 3536]). 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.,
(21) 
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
(22) 
i.e., the approximation converges to in the meansquares sense (cameron_1947, ). Basis functions for many common PDFs are well known as part of the WienerAskey 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
(23) 
with the th trigonometric moments ()
(24)  
(25) 
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,
(26)  
(27) 
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
(28) 
and the circular standard deviation is
(29) 
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. 5051) and the von Mises distribution (vonmises_1918, ). The (onedimensional) Wrapped Normal Density (WND)
(30) 
is parameterized by the mean direction and the variancelike 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
(31) 
Initially proposed for applications in atomic physics, the von Mises distribution (VMD) is
(32) 
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
(33) 
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,
(34)  
(35)  
(36) 
where
(37) 
and are dubbed the Verblunsky coefficients. These coefficients are unique for a given OPUC defined by measure in Eq. (21), and
(38) 
For the form of the characteristic function in Eq. (23), the Toeplitz matrix for a given measure on is
(39) 
As described in (szego_1975, , pp. 287288),
(40) 
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 finiteprecision 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 RogersSzegő Polynomials
Originally presented by Gabor Szegő (szego_1926, ) and based on the Hermite polynomials of Leonard Rogers (rogers_1893, ; rogers_1894, ), the RogersSzegő polynomials are orthogonal with respect to . The (normalized) RogersSzegő polynomial of degree is (atakishiyev_1994b, )
(41) 
with parameter and the binomial coefficients
(42) 
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 RogersSzegő polynomials are
(43) 
These polynomials form an orthogonal basis with regards to the WND measure (e.g., see discussion in simon_2004 (); atakishiyev_1994b ()) given input argument
(44) 
and parameter
(45) 
Note that, via Euler’s formula, the RogersSzegő polynomials are trigonometric polynomials of the random angle . Figure 2 illustrates the real and imaginary components of the RogersSzegő polynomials for , , and . Note that the polynomials are periodic in .
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
(46)  
(47) 
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 RogerSzegő 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 ,
(48)  
(49)  
(50)  
(51) 
since and for . To solve for the term of the characteristic function,
(52)  
(53)  
(54) 
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.
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
(57)  
(58)  
(59) 
and the orthogonality of leaves the function solely dependent on the zeroth order coefficients and the normalizing constants. Similarly, when considering the second moment,
(60)  
(61)  
(62) 
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 orbitstate 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
(63) 
with random decayrate coefficient and deterministic solution
(64) 
For this case, , where is a random variable on the unit circle, and time unit. The mean and variance of are
(65)  
(66) 
which, for a given density of , are computed numerically via quadrature using 1000 uniformly spaced nodes on the unit circle^{1}^{1}1See CIRCLE_RULE software found at https://people.sc.fsu.edu/jburkardt/c_src/circle_rule/circle_rule.html [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 surrogatederived statistics are compared to numeric approximations of Eqs. (65) and (66). Let and denote the PCE coefficients for and , respectively. In the following, multiindices 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
(67) 
where
(68) 
The initial conditions are
(69) 
and the PCE for has coefficients
(70) 
which is the analytic solution for a PCE such that . The ODE for the PCE coefficients (see Eq. (67)) may be used in any RungeKutta 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 precomputed using quadrature integration with the appropriate polynomials and stored for use in propagating the PCE coefficients.
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
(71) 
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 nontrivial 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  
VMD  
Concentrated  WND  
VMD 
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 nonzero 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 RogersSzegő polynomials, do not suffer from this issue.
v.2 Earth Orbit Cases
Of particular interest in this work is propagation of orbitstate 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 twobody 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 DormandPrince 8(7) prince_1981 () embedded RungeKutta propagator with a relative tolerance of . The gravitational force model uses km/s and . For the sake of simplicity, transformation between the Earthfixed 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 Earthcentered orbit propagation for the sake of simplicity and duplication, previous demonstrations of surrogate methods for uncertainty propagation consider higherfidelity dynamics jones_2013a (); jones_2015b (); balducci_2017 ().
Moment  

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 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 nonGaussian 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 RootMeanSquare (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 twobody only dynamics to propagate the angle:
(72) 
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 RogerSzegő 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.
Method  Basis  QOI  Rel. Error Mean  Rel. Error STD  RMS Error [deg] 

PCE  Hermite  
Hermite  
RogersSzegő  
SR  Hermite  
RogersSzegő 
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 RogersSzegő) with QOI exhibit only slight differences in the RMS error for the PCE surrogates. Otherwise, no significant difference is exhibited.
Figure 6 illustrates the effect that QOI selection has on the coefficients of the PCE. When using RogersSzegő 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 RogersSzegő results. Using the RogerSzegő 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.
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: SMAOnly 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. 653654) 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 RogersSzegő polynomials, and are compared to a Monte Carlo analysis with samples and a resulting precise to approximately four digits. Figure 8 presents the posterior (nonGaussian) marginal PDF for using these Monte Carlo samples.
Figure 9 presents the accuracy of the PCEdetermined 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.
Figure 10 demonstrates the RMS error as a function of . Since the PCEproduced 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 crossvalidation 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.
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 nonGaussian, which results from shortperiod 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 RogersSzegő 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.
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 nonGaussian. 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.
QOI ()  WND  WND  VMD 

deg  deg  
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 polynomialbased methods require a basis orthogonal with respect to the input probability density functions. Random inputs described by the wrapped normal density employ the RogersSzegő 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 orbitstate uncertainty when using the equinoctial orbital elements.
Acknowledgements.
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.References
References
 (1) Junkins, J. L., Akella, M. R., and Alfriend, K. T., “NonGaussian 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 TwoBody Problem,” Journal
of Guidance, Control, and Dynamics, Vol. 35, No. 2, 2012, pp. 497–509.
doi:10.2514/1.54385. 
(3)
DeMars, K. J., Bishop, R. H., and Jah, M. K., ‘‘EntropyBased Approach
for Uncertainty Propagation of Nonlinear Dynamical Systems,” Journal of
Guidance, Control, and Dynamics, Vol. 36, No. 4, 2013, pp. 1047–1057.
doi:10.2514/1.58987.  (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 NonIntrusive Polynomial Chaos,” Journal of
Guidance, Control, and Dynamics, Vol. 36, No. 2, 2013, pp. 430–444.
doi:10.2514/1.57599. 
(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. 12, 2017,
pp. 105–136.
doi:10.1007/s1056901797677. 
(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.
doi:10.2514/1.58068. 
(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.
doi:10.1016/j.asr.2013.08.027. 
(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.
doi:10.2514/1.G000595.  (10) Feldhacker, J. D., Smith, J., Jones, B. A., and Doostan, A., “MultiElement Trajectory Models for Satellite Tour Missions,” in “AIAA/AAS Astrodynamics Specialist Conference,” AIAA 20165502, 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.
doi:10.1137/130917296.  (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 RedHorse, J., “Propagation of probabilistic
uncertainty in complex physical systems using a stochastic finite element
approach,” Physica D: Nonlinear Phenomena, Vol. 133, No. 14, 1999,
pp. 137–144.
doi:10.1016/S01672789(99)001025. 
(17)
Ghanem, R. G., “Ingredients for a General Purpose Stochastic Finite
Elements Implementation,” Computer Methods in Applied Mechanics and
Engineering, Vol. 168, No. 14, 1999, pp. 19–34.
doi:10.1016/S00457825(98)001066. 
(18)
Xiu, D. and Karniadakis, G. E., “The WienerAskey Polynomial Chaos
for Stochastic Differential Equations,” SIAM Journal of Scientific
Computing, Vol. 24, No. 2, 2002, pp. 619–644.
doi:10.1137/S1064827501387826.  (19) Doostan, A., Iaccarino, G., and Etemadi, N., “A leastsquares approximation of highdimensional uncertain systems,” Tech. Rep. Annual Research Brief, Center for Turbulence Research, Stanford University, 2007.

(20)
Doostan, A. and Iaccarino, G., “A leastsquares approximation of
partial differential equations with highdimensional random inputs,”
Journal of Computational Physics, Vol. 228, No. 12, 2009, pp.
4332–4345.
doi:10.1016/j.jcp.2009.03.006. 
(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.
doi:10.1016/j.jcp.2009.01.006. 
(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.
doi:10.1007/s1183101090541. 
(23)
Doostan, A. and Owhadi, H., “A Nonadapted Sparse Approximation of
PDEs with Stochastic Inputs,” Journal of Computational Physics,
Vol. 230, No. 8, 2011, pp. 3015–3034.
doi:10.1016/j.jcp.2011.01.002. 
(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.
doi:10.1016/j.jcp.2011.10.028. 
(25)
Doostan, A., Validi, A., and Iaccarino, G., “Nonintrusive lowrank
separated approximation of highdimensional stochastic models,”
Computation Methods in Applied Mechanical Engineering, Vol. 263, 2013,
pp. 42–55.
doi:10.1016/j.cma.2013.04.003. 
(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.
doi:10.1137/110854059. 
(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.
doi:10.1002/nme.4761. 
(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.
doi:10.1137/130949154. 
(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.
doi:10.1137/130929461. 
(30)
Gerritsma, M., van der Steen, J.B., Vos, P., and Karniadakis, G.,
“Timedependent generalized polynomial chaos,” Journal of
Computational Physics, Vol. 229, No. 22, 2010, pp. 8333 – 8363.
doi:10.1016/j.jcp.2010.07.020.  (31) Wan, X. and Karniadakis, G. E., “An adaptive multielement 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., “MultiElement Generalized Polynomial
Chaos for Arbitrary Probability Measures,” SIAM Journal on Scientific
Computing, Vol. 28, No. 3, 2006, pp. 901–928.
doi:10.1137/050627630.  (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 NonIntrusive Polynomial Chaos Method for Uncertainty Propagation in CFD Simulations,” in “44th AIAA Aerospace Sciences Meeting and Exhibit,” AIAA 2006891, 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.
doi:10.1137/070710524.  (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
NonLinear Functionals in Series of FourierHermite Functionals,”
Annals of Mathematics, Vol. 48, No. 2, 1947, pp. 385–392.
doi:10.2307/1969178. 
(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.
doi:10.1098/rspa.1953.0064.  (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 FisherBingham 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.
doi:10.1080/02331880500031597. 
(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.
doi:10.1111/j.1467842X.1981.tb00763.x.  (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.
doi:10.1137/0903018.  (49) Szegő, G., “Ein Beitrag zur Theorie der Thetafunktionen,” Sitzungsberichte der Preussischen Akademie der Wissenschaften, physikalischmathematische, pp. 242–252.

(50)
Rogers, L. J., “Second Memoir on the Expansion of certain Infinite
Products,” Proceedings of the London Mathematical Society, Vol. s125,
No. 1, 1893, pp. 318–343.
doi:10.1112/plms/s125.1.318. 
(51)
Rogers, L. J., “Third Memoir on the Expansion of certain Infinite
Products,” Proceedings of the London Mathematical Society, Vol. s126,
No. 1, 1894, pp. 15–32.
doi:10.1112/plms/s126.1.15. 
(52)
Atakishiyev, N. M. and Nagiyev, S. M., “On the RogersSzegő
polynomials,” Journal of Physics A: Mathematical and General, Vol. 27,
No. 17, 1994, pp. L611–L615.
doi:10.1088/03054470/27/17/003. 
(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.
doi:10.2307/2346732.  (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 RungeKutta
formulae,” Journal of Computational and Applied Mathematics, Vol. 7,
No. 1, 1981, pp. 67–75.
doi:10.1016/0771050X(81)900103.  (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.