Satellite collision probability estimation using polynomial chaos expansions

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research 52 (2013) 1860–1875 www.elsevier.com/locate/asr

Satellite collision probability estimation using polynomial chaos expansions Brandon A. Jones ⇑, Alireza Doostan Department of Aerospace Engineering Sciences, University of Colorado Boulder, UCB 431, Boulder, CO, United States Received 10 April 2013; received in revised form 16 August 2013; accepted 23 August 2013 Available online 4 September 2013

Abstract This paper presents the application of polynomial chaos (PC) to estimating the probability of collision between two spacecraft. Common methods of quantifying this probability for conjunction assessment use either Monte Carlo analyses or include simplifying assumptions to improve tractability. A PC expansion, or PCE, provides a means for approximating the solution to a large set of stochastic ordinary differential equations, which includes orbit propagation. When compared to Monte Carlo methods, non-intrusive, i.e., sampling-based, PCE generation techniques may greatly reduce the number of orbit propagations required to approximate the possibly non-Gaussian, a posteriori probability density function. The presented PC-based method of computing collision probability requires no fundamental simplifying assumptions, and reduces the computation time compared to Monte Carlo. This paper considers two cases where the common conjunction assessment assumptions are no longer valid. The results indeed demonstrate a reduction in computation time when compared to Monte Carlo, and improved accuracy when compared to simplified techniques. Ó 2013 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Conjunction assessment; Uncertainty propagation; Polynomial chaos; Space situational awareness; Orbital debris

1. Introduction As the space in low-Earth and geosynchronous orbit becomes increasingly congested, rigorous methods of estimating the risk of collision between any two objects becomes increasingly important. Such methods must be highly accurate to: (i) identify possible collisions, and (ii) minimize the number of false alarms that waste fuel and disrupt satellite operations. This paper considers a new application of polynomial chaos (PC) methods (Ghanem and Spanos, 1991; Ghanem and Dham, 1998; Ghanem, 1999; Ghanem and Red-Horse, 1999; Xiu and Karniadakis, 2002; Le Maitre and Knio, 2010) for the quantification of collision risks in scenarios where established methods no longer apply. Several methods of computing a collision probability already exist and are used for satellite operations today. ⇑ Corresponding author. Tel.: +1 303 492 3753.

E-mail addresses: [email protected] (B.A. Jones), Alireza. [email protected] (A. Doostan).

The most commonly discussed methods (Alfriend et al., 1999; Patera, 2001; Alfano, 2005) share several assumptions: A1. the probability density function (pdf) describing the position state is a multivariate Gaussian distribution, A2. the change in relative position between the two objects is linear with constant velocity over the period of conjunction, A3. uncertainties in the velocities may be ignored, A4. the uncertainties for the two objects are uncorrelated, A5. the position uncertainties remain unchanged over the period of conjunction, and A6. both spacecraft are described as spherical bodies. This paper considers situations where at least one of these conditions is not satisfied and these classic methods fail to accurately approximate the probability of collision. Other methods of collision probability estimation have been developed to mitigate one or more of the common assumptions, but not all. For example, Patera (2003), Chan

0273-1177/$36.00 Ó 2013 COSPAR. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.asr.2013.08.027

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

(2004) and Alfano (2006) discuss methods to compute the probability of collision in the presence of nonlinear motion, but do not consider non-Gaussian state uncertainties. Recent research demonstrates that the Gaussian pdf assumption is not strictly true in astrodynamics, with various methods suggested for non-Gaussian uncertainty propagation (Jones et al., 2013; Fujimoto et al., 2012; Fujimoto and Scheeres, 2012; DeMars et al., 2009; Horwood et al., 2011; Sabol et al., 2010; Park and Scheeres, 1367; Valli et al., 2013). The method based on PC discussed in this paper makes no assumption on the a posteriori state pdf, and includes a test case where the position states for both satellites are non-Gaussian. The increasing interest in formation flying (e.g, see Sharma and Curtis, 2005) and fractionated spacecraft,1 offers a new challenge in the form of conjunction assessment between satellites in a single constellation. Navigation systems for such formations may use relative measurements, e.g., cross-link range, to improve relative-state accuracy, which introduces correlations between satellite state pdfs (e.g., Carpenter, 2002). The uncorrelated uncertainty assumption (A4) is introduced in the classic collision probability estimation methods when forming the relative state pdf (mean and covariance) from absolute state information. While using a direct relative state navigation solution bypasses this assumption and is the ideal input for computing the collision probability, such quantities of interest must be propagated to the time of closest approach with sufficient planning time to allow for collision avoidance maneuver design. For nonlinear relative dynamics, the relative position pdf may be non-Gaussian during the conjunction. This violates two of the previously discussed assumptions (A1 and A2) for the classic collision probability estimation methods. Hence, the classic methods may not provide sufficient accuracy for potential collisions between satellites in such constellations. This paper presents a method of propagating the correlated absolute state uncertainties to compute collision probabilities for formation-flight cases. Conjunction assessment, specifically collision probability estimation, based on properly formulated Monte Carlo (MC) methods provide the most accurate, but computationally expensive, estimate of the collision risk. Such methods remove the previously stated assumptions, and are only limited by the accuracy of the propagator, the realism of the initial uncertainty, and computation time. Additionally, a recent assessment of the astrodynamics standards employed by the United States Air Force Space Command acknowledges that particle-based methods are required for some cases (National Research Council, 2012), motivating the development of such techniques. However, the computation time requirements can pose a major limitation in adoption of such methods (Sabol 1

For example, see a description of the United States Defense Advance Research Projects Agency (DARPA) F6 program at http://www.darpa.mil/Our_Work/tto/Programs/System_F6.aspx [accessed March 2013].

1861

et al., 2011; Dolado et al., 2011; Garmier et al., 2011). Several tools exist to reduce the long execution times associated with Monte Carlo methods, with these tools leveraging off of both improvements in computation capabilities and theories developed to reduce the number of MC trials. Sabol et al. (2011) demonstrate the use of a parallelized algorithm to perform Monte Carlo analyses within seconds per satellite. Dolado et al. (2011) consider importance sampling to reduce the number of orbit propagations required to quantify collision probability. Garmier et al., 2011 reduce the overall computation time by filtering out scenarios with a sufficiently low probability of collision, thereby limiting the total number of Monte Carlo samples for a given space object catalog. The PC-based methods discussed in this paper share similarities with those based on Monte Carlo, but use polynomial surrogates to reduce the total number of orbit propagations. Use of these polynomial surrogates makes the generation of realizations at a future point in time more tractable, with thousands of evaluations per second possible on a single computer processor. Parallelized implementations of the PC methods provide further reductions in computation time. A polynomial chaos expansion (PCE) allows for the approximation of the solution to a stochastic differential equation that is square-measurable, possibly non-Gaussian, with respect to the input uncertainties (Ghanem and Red-Horse, 1999; Ghanem and Spanos, 1991; Ghanem and Dham, 1998; Ghanem, 1999; Xiu and Karniadakis, 2002, 2003; Le Maitre et al., 2004; Babusˇka et al., 2005; Wan and Karniadakis, 2005; Xiu, 2009, 2010; Le Maitre and Knio, 2010). PC-based methods use a projection of the stochastic solution at a given time tf onto a basis of orthogonal polynomials in stochastic variables that is dense in the space of finite-variance random variables at some initial time t0 . After generating a PCE, the polynomial representation provides a response surface in the space of the random inputs, thereby allowing for the generation of realizations at time tf when given random inputs at the initial time. In generating these realizations, no additional evaluations of the orbit propagator are required, thereby reducing the total computation time. The current authors already presented the use of PC for propagating orbit uncertainty over relatively long durations (10 days), and results demonstrate the ability to accurately approximate a non-Gaussian pdf at the final time (Jones et al., 2013). Additionally, Jia et al. (2012) briefly demonstrated the use of PC in the two-body problem. The current paper is a follow-up to a previous conference paper that only considered an instantaneous probability of collision (Jones et al., 2012), and extends the methods to quantifying a satellite survival probability. This research considers the use of PC, specifically the generated response surfaces, to estimate collision probabilities greater than or equal to 105 . In conjunction assessment operations for manned missions, the United States National Aeronautics and Space Administration (NASA) uses a decision threshold of 105 for avoidance maneuver

1862

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

execution, and increases the threshold to 104 for robotic missions.2 As discussed later, the presented PC-based methods do not require the common limiting assumptions to yield an accurate estimate of collision probability. However, the current paper does not investigate several important elements of the larger conjunction assessment problem. In particular, the presented methods assume that a potential collision has already been identified. The work of Hoots et al. (1984), Woodburn and Dichmann (1997), Woodburn et al. (2009), Alfano (2012), Armellin et al. (2012) and Mercurio et al. (2012) have developed methods for identifying a conjunction. Additionally, the current paper assumes some time period of the conjunction has been provided. See Alfano and Negron (1993), Alfano (1994) and Alfano and Greer (2003) for possible techniques. PC-based methods for these tasks will be investigated in future work. The rest of this paper is organized as follows. Section 2 provides an overview of PC and the associated methods employed in this paper. Section 3 then discusses how methods based on PC and Monte Carlo may be used to estimate the probability of collision. Two test cases are then considered in Section 4 to demonstrate the computation savings when using PC. Finally, Section 5 provides conclusions and summarizes future work identified in the paper.

tions (ODEs) that, for the present case, describe the temporal evolution of a satellite’s state as Aðt; n; Xðt0 ÞÞ ¼ 0;

ðt; nÞ 2 ½t0 ; tf   Cd ;

ð1Þ where A is a stochastic ordinary differential operator, t 2 ½t0 ; tf  is the temporal variable, and Xðt0 Þ is the initial condition at t0 . This work considers the propagation of orbit state uncertainty using PC with the goal of estimating collision probability. Upon defining the uncertainty in Xðt0 Þ 2 Rn on the probability space ðS; F ; PÞ, PC-based methods are then used to generate a solution to the stochastic operator A that describes the space of possible solutions Xðtf Þ at any future time tf . We note that n and d, i.e., the dimension of the approximated state vector and the stochastic dimension, are not necessarily equal. PC methods use the spectral basis fwa g, which is orthonormal with respect to the probability density qðnÞ of the random inputs n, i.e., Z wa ðnÞwb ðnÞqðnÞ dn ¼ dab ; ð2Þ Cd

where dab is the Kronecker delta function, and a and b are multi-indices. Given a multi-index a ¼ ða1 ; a2 ; . . . ; ad Þ, wa ðnÞ ¼ wa1 ðn1 Þwa2 ðn2 Þ . . . wad ðnd Þ

2. Polynomial chaos Methods based on polynomial chaos provide a means for generating an approximation to the solution of a stochastic system by projecting it onto a basis of spectral polynomials of random inputs. Wiener (1938) first proposed this type of approximation, with methods based on Hermite chaos expansions first established (Ghanem and Spanos, 1991; Ghanem and Dham, 1998; Ghanem, 1999; Ghanem and Red-Horse, 1999) and later generalized to non-Gaussian input uncertainties (Xiu and Karniadakis, 2002; Soize and Ghanem, 2004). The generated PCE represents the stochastic solution as a linear expansion of multivariate polynomials that are functions of the input uncertain variables n 2 Rd where d is the stochastic dimension of the system. Solving for the PCE means estimating the coefficients, or weights, of the expansion. PC-based methods provide a computationally efficient means to represent any square-measurable, possibly non-Gaussian distribution, and have already been demonstrated for orbit propagation (Jones et al., 2013). Let ðS; F ; PÞ be a probability space with the sample space S and the probability measure P on r-algebra F . The d random inputs to the system, assumed here to be independent, are n 2 Rd : S ! Cd # Rd on ðS; F ; PÞ, where these elements characterize the input uncertainties. Additionally, we denote the set of ordinary differential equa2 See page 17 of eBook at http://www.nasa.gov/offices/oce/appel/ knowledge/publications/appel-releases-ibook.html [accessed March 2013].

P  a:s: in S;

ð3Þ

is the tensor product of univariate polynomials wai where the non-negative integer ai specifies the degree of wai ðnÞ. In other words, the multivariate polynomials wa ðnÞ are simply the product of univariate polynomials wai ðni Þ. Given n and qðnÞ, Eq. (2) specifies the basis fwa g used in the PCE. The current paper only considers n 2 N ð0; I d Þ, i.e., the random inputs are multivariate Gaussian with zero mean and a d  d identity matrix I d for the variance-covariance matrix. For this selection of n and qðnÞ; wa are Hermite polynomials (Ghanem and Spanos, 1991). Although n are Gaussian distributed, the uncertainty of Xðt0 Þ is not necessarily Gaussian. Additionally, the elements of n do not have to be identically distributed, i.e., random inputs may combine, for example, a mixture of Gaussian and uniform distributions (Soize and Ghanem, 2004). In the context of PC, the orbit solution Xðt; nÞ to Eq. (1) is then approximated by the finite series X ^ nÞ ¼ ca ðtÞwa ðnÞ; ð4Þ Xðt; a2Kp;d

^ signifies the approximate solution to the where the hat in X true value X, and the vector of coefficients ca ðtÞ 2 Rn yields the PCE solution at t. The finite set Kp;d :¼ fða1 ; . . . ; ad Þ : kak1 6 p; kak0 6 dg

ð5Þ

is the set of multi-indices of size d defined on the nonnegative integers and yielding a polynomial basis of total degree not larger than p. The PCE approximation is then refined by increasing p to achieve a target accuracy. If, for any t 2 ½t0 ; tf ; Xðt; nÞ 2 L2 ðCd Þ, then the PCE

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

approximation of Eq. (4) converges in the L2 sense as p ! 1 (e.g., Xiu, 2010). We note that a vector representation for ca effectively allows for the simultaneous representation of n independent PCE’s, i.e., one for each element of the vector X 2 Rn . The number of terms in Eq. (4) is P :¼

ðp þ dÞ! : p!d!

ð6Þ

As implied in this equation, P increases exponentially with p and d, leading to the issue of curse-of-dimensionality. While this can impose computational difficulties when d or p is large, ongoing research seeks to reduce the dimensionality and computation costs for such cases (Doostan et al., 2007, 2013; Doostan and Iaccarino, 2009; Doostan and Owhadi, 2011). The PCE in Eq. (4) provides a method for generating an ^ f ; nÞ when given a random approximate realization Xðt input n, i.e., it describes the response of Xðtf Þ to the input n. In this way, the PCE yields a solution to a stochastic system, as opposed to propagating the state pdf. Information on the a posteriori pdf may be inferred from the stochastic solution via sampling of the polynomial surrogate or through analytic functions of ca (e.g., Ghanem and Spanos, 1991; Xiu, 2010). Given sufficient accuracy of ca , the polynomial surrogates also provide a means for bypassing an ODE solver and approximating the Monte Carlo results using computationally inexpensive polynomial evaluations. This latter property allows for the use of PC in conjunction assessment, which is described further in Section 3. The remainder of this section describes a computationally efficient method for generating a PCE in the context of astrodynamics. 2.1. PCE solution methods This study only considers a regression-based method to generate a PCE approximating the stochastic solution of the propagated orbit at time t. Such methods treat already existing ODE solvers as a black box, and require no alterations to existing orbit propagation software, hence they are dubbed non-intrusive. This allows us to wrap existing software to propagate the stochastic solution and enables their use in conjunction assessment with minimal alteration to existing tools. The algorithm behind the non-intrusive methods may be summarized by:

Several things should be noted about this procedure. First, PCE solutions at intermediate steps, e.g., integration time steps, may be easily generated using the same samples at those times. More samples are only required if p increases at the intermediate time tk . Hence, ca ðtÞ may be generated at multiple times using the same samples. Second, there are several methods of computing ca ðtf Þ in step 4 above using realizations Xðtf ; ni Þ, including, but not limited to: least squares regression (Hosder et al., 2006), pseudo-spectral collocation (Xiu, 2010; Le Maitre and Knio, 2010), Monte Carlo sampling (Le Maitre and Knio, 2010), and compressive sampling (Doostan and Owhadi, 2011). Each of these methods have several advantages and disadvantages which will not be discussed here. The current study considers solutions based on regression, which were previously demonstrated in Jones et al. (2013). The next section provides a brief overview of the regression methods. 2.1.1. PCE solutions via regression The method to solve for ca based on regression uses random samples ni from the density function qðnÞ. Given the M PC propagated states, one solution for the coefficients ca minimizes the cost function PC 1X T i ; 2 i¼1 i

M

J ðca Þ ¼

ð7Þ

where ^ f ; ni Þ  Xðtf ; ni Þ: i ¼ Xðt

ð8Þ

Eq. (7) is simply the cost function of the least squares estimator. Upon inspection of Eq. (4), we may represent the PC PCE evaluation of the M PC random inputs fni gM i¼1 as the linear system 32 T 3 2 T 3 2 ca1 ðtf Þ wa1 ðn1 Þ    waP ðn1 Þ X ðtf ; n1 Þ 76 7 6 7 6 6 wa1 ðn2 Þ    waP ðn2 Þ 76 cTa2 ðtf Þ 7 6 X T ðtf ; n2 Þ 7 76 7 6 7 6 76 . 7 ¼ 6 7 6 .. .. .. .. 7 7 7 6 6 6 . . . . . 54 . 5 4 5 4 T T wa1 ðnM PC Þ    waP ðnM PC Þ X ðtf ; nM PC Þ caP ðtf Þ ð9Þ where the right-hand-side is populated by the propagated samples Xðtf ; nÞ. Alternatively, one may write Eq. (9) as HC tf ¼ Y tf ;

1. Generate M PC realizations ni . For the current paper, these samples are based on a random sampling of the distribution N ð0; I d Þ. 2. For each of the M PC random vectors, use the initial uncertainty in Xðt0 ; ni Þ to generate a realization based on the random input ni . 3. Propagate Xðt0 ; ni Þ to tf using an existing ODE solver to obtain Xðtf ; ni Þ for each of the M PC realizations. 4. Solve for the PC coefficients ca ðtf Þ in Eq. (4) based on Xðtf ; ni Þ and the method of choice. This yields ^ f ; nÞ. a PCE representation Xðt

1863

ð10Þ

where H is the M PC  P measurement matrix on the left hand side of Eq. (9), C tf is the matrix of PCE coefficients used to approximate Xðtf ; ni Þ, and Y tf is comprised of the propagated states at tf . Given this formulation, ^ t ¼ ðH T HÞ1 H T Y t ; C f f

ð11Þ

which is the solution to the normal equation associated with the least squares estimator. This method of solving for the PCE coefficients requires M PC P P measurements. Also, H does not change as a function of tf . This allows

1864

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

for a computationally efficient method for generating PCE coefficients at multiple points in time, which is described in Section 2.1.3.

extend the process to a K-fold cross validation approach. However, this algorithm proves sufficient for the tests described in Section 4.

2.1.2. Generation of the PCE approximation The previous section described the algorithm to generate a PCE using regression. However, this method requires the selection of M PC and p. Characterizing the approximation accuracy of a non-intrusive PCE solution requires analysis of the estimation residuals based on a comparison with ^ t . Since samples independent of those used to compute C k this work focuses on PCE’s used for conjunction assessment, validation of the approximation focuses on the root-mean-square (RMS) residuals of the Cartesian position rðtÞ. Although variations on this technique exist, the method considered in this work iteratively increases M PC for an assumed p until the regression-based approximation achieves a specified convergence. To begin the PCE generM0 ation process, the set of samples fni ; Xðtk ; ni Þgi¼1PC are generated where

2.1.3. PCE generation and evaluation at multiple times The PC-based methods of conjunction assessment require the generation of realizations at some time tk . In some cases, analysis may be required at multiple tk , where k ¼ 1; . . . ; K. This section describes computationally efficient methods for generating these solutions. Solving for a PCE at multiple times tk using regression may be performed more efficiently using linear algebra operations. As indicated in Eq. (9), H is only a function of the number of inputs M PC ; P , and the random vectors ni . Assuming these inputs are fixed, then the same H may be employed to generate PCE’s at all times. Specifically, given the propagated solution at all times tk ,     ^t  C ^ t ¼ H T H 1 H T ½ Y t1    Y tk  ð17Þ C K 1

M 0PC ¼ M PC þ M V ; M PC ¼ ceilðcP Þ;     P M V ¼ max ceil ; 10 ; 10

ð12Þ ð13Þ ð14Þ

where c is a design parameter, ceilðÞ denotes rounding up to the next integer, and M V is the number of validation samples used in the procedure. In this paper, c ¼ 1:1, and Eq. (14) ensures the consideration of at least ten validation samples. Selection of all M 0PC samples is based on random draws from N ð0; I d Þ, and the first ceilðcP Þ samples are ^ t , while the M V validation samples are used to compute C k saved to test for solution convergence. Using the validation sample inputs ni to generate the PCE-determined position vector ^rðtk ; ni Þ, the position RMS residual sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PM V i¼1 Dri  Dri ; ð15Þ RMS ¼ MV where Dri ¼ ^rðtk ; ni Þ  rðtk ; ni Þ

ð16Þ

is computed and compared to a user-specified convergence tolerance. If RMS fails to meet the specified convergence criteria, more samples are generated (increasing the total ^ t and RMS . number M 0PC ) and used to compute a new C k This process continues until the convergence criteria are satisfied, which yields a PCE with a position RMS accuracy of approximately RMS . The PCE accuracy resulting from this procedure is verified in Section 4. The value of p is a design parameter for the PCE generation process and fixed at process initialization. Since M V additional samples are required for solution validation, the number M 0PC denotes the total number of random trials required to generate the PCE. Follow-up work will extend the iterative solver described previously to include the selection of p and

yields the PCE coefficients for all times of interest. Eq. (4) may also be reformulated to allow for M evaluations of the PCE at multiple times when given the set of M random variables fni gi¼1 . Assuming the coefficient matrix ^ ^ k ; ni Þ based on C tk is already known, then realizations Xðt M fni gi¼1 may be generated via the linear equation 3 2 ^ T ðtk ; n1 Þ X 7 6 ^T 6 X ðtk ; n2 Þ 7 7 6 ^t : ^ ð18Þ Y tk ¼ 6 7 ¼ HC k .. 7 6 . 5 4 ^ T ðtk ; nM Þ X ^ t at multiple times When given PCE coefficient matrices C k t1 ; t2 ; . . . ; tK , then     ^t  C ^t : ^t  Y ^t ¼ H C ð19Þ Y K K 1 1 Already existing tools provide parallel implementations of Eqs. (17)–(19). For example, some variations of the Basic Linear Algebra Subprograms (BLAS)3 tools, e.g., PBLAS,4 already support multiprocessor computing systems. Hence, implementation of PC on, for example, a graphics processing unit (GPU)5 requires minimal software development. 2.2. PCE coordinate system selection considerations As demonstrated in the literature, orbit elements provide a better coordinate system for propagating orbit uncertainty since the pdf remains Gaussian for longer time periods (Jones et al., 2013; Sabol et al., 2010; Fujimoto and Scheeres, 2012; Fujimoto et al., 2012). As demonstrated in Jones et al., 2013, using an element set in the generation of 3

BLAS Website: http://www.netlib.org/blas/ [Accessed March 2013]. PBLAS Website: http://www.netlib.org/scalapack/pblas_qref.html [Accessed March 2013]. 5 BLAS for GPUs: https://developer.nvidia.com/cublas [Accessed March 2013]. 4

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

the PCE reduces d, and, thus, M PC P P . However, the numerical results in Section 4 use a Cartesian coordinate system. These tests are based on cases provided in the literature, which are defined in Cartesian coordinates. The initial covariance matrix may be mapped into an element set, e.g., the equinoctial elements, using the unscented transformation, but at the sacrifice of pdf accuracy due to the nonlinear transformation. Unfortunately, using a pdf representation of an orbit element set poses a difficulty in performing conjunction assessment: the mapping from elements to Cartesian position is not injective. Specifically, two distinct element vectors may describe a satellite state at the same position but with different velocities. Hence, when performing conjunction assessment, any description of the state pdf must be converted to Cartesian coordinates before computing the probability of collision. Such methods may yield a loss of accuracy in the nonlinear transformation of the pdf. In the case of PC-based methods, equinoctial element realizations generated using Eq. (4) or (19) may be converted to Cartesian coordinates for calculation of the empirical probability of collision. Transforming such PCE realizations for use in conjunction assessment removes the inaccuracies introduced in the element-to-Cartesian pdf transformation. Although such methods are not demonstrated in this paper, they still apply with this extra transformation.

The methods employed for estimating the probability of collision using PC are similar to Monte Carlo techniques, but rely on the polynomial surrogates of Eq. (4) to reduce the overall computation time. The following discussion begins with a description of Monte Carlo methods, which then serves as a framework for presenting the PC-based techniques. The presentation also assumes an initial Gaussian distribution, but may be altered for any arbitrary pdf at t0 . 3.1. Monte Carlo estimation of P c For the state vector  T X l ðt0 Þ ¼ rTl ðt0 Þ r_ Tl ðt0 Þ

ð20Þ

comprised of the position rl and velocity r_ l for satellite l, let  l ðt0 Þ be the mean at time t0 with covariance matrix Pl . In X a Monte Carlo test, M MC trials are generated using, for example,  l ðt0 Þ þ Ll ni X l ðt0 ; ni Þ ¼ X

satellites included in the analysis. Given two satellites (l = 1, 2), their separation at time tk is sðtk ; ni ; n0i Þ ¼ kr1 ðtk ; ni Þ  r2 ðtk ; n0i Þk2

ð21Þ

where Ll is the lower-triangular Cholesky decomposition of Pl ; i ¼ 1; . . . ; M MC , and ni  N ð0; I 6 Þ. Each trial X l ðt0 ; ni Þ is then propagated to some time tk to yield X l ðtk ; ni Þ. This procedure is conducted for any number of

ð22Þ

where the prime on n0i is only meant to indicate a different set of samples, i.e., ni –n0i , and krk2 denotes the L2 norm of r. The empirical, instantaneous probability of collision is then P c ðRÞ ¼

countðsðtk ; ni ; n0i Þ 6 RÞ ; M MC

ð23Þ

where R is a given collision hard-body radius, and the count() operator indicates the number of true results of the argument over i ¼ 1; . . . ; M MC . Except when considering a specific value of R, the instantaneous probability of collision is denoted by P c . Using the propagated realizations at multiple points in time, i.e., Xðtk ; ni Þ for tk ; k ¼ 1; . . . ; K, a probability of survival P s may be computed. The complement of the survival probability, P 0s , describes the probability of collision over a period of time as opposed to a single epoch. To compute P s or P 0s , let mðni ; n0i Þ ¼ min sðtk ; ni ; n0i Þ: k¼1;...;K

ð24Þ

This yields the minimum distance between the two satellites over the time period ½t1 ; tK . Hence, P 0s ðRÞ ¼ 1  P s ðRÞ ¼

3. Collision probability estimation

1865

countðmðni ; n0i Þ 6 RÞ : M MC

ð25Þ

Some use the term cumulative probability to describe P 0s (e.g., Alfano, 2009), but such terminology can imply a sum of instantaneous probabilities and is avoided in this paper. The method defined in Eqs. (24) and (25) assumes the times tk are dense enough to provide adequate accuracy when computing P 0s . Alternatively, interpolation may be employed with a given minimization algorithm to improve the precision of the time of minimum distance, and, thus, increase accuracy. This would come at the cost of additional computation time and is not used in the current paper. The presentation of Eqs. (20)–(25) assumes the position vectors for the two satellites are uncorrelated. However, only a slight alteration to the algorithm is required to account for cross-satellite correlations. For the combined state vector of satellites 1 and 2,  T ð26Þ X 1;2 ¼ rT1 ðt0 Þ r_ T1 ðt0 Þ rT2 ðt0 Þ r_ T2 ðt0 Þ ; let N ðX 1;2 ; P1;2 Þ describe the joint pdf of the two satellite states. Then a single input vector ni 2 R12 is required to generate the realizations at t0 . In this case, ni ¼ n0i in Eqs. (22)–(25), i.e., the realizations are a function of the same random inputs. Otherwise, the algorithm is the same. For two satellites with uncorrelated states, one may perform a one-to-one or all-on-all comparison of the Monte Carlo realizations. The one-to-one methods use M MC realizations for each satellite, with each one employed only

1866

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

once in Eqs. (22) and (23) (e.g., Ghrist and Plakalovic, 2012). Alternatively, an all-on-all comparison considers M 2MC possible combinations, implying that fewer propagated samples are required to yield a given number of comparisons (e.g., Sabol et al., 2011). In this latter case, not every comparison is independent, hence it qualifies as an approximation method to reduce the number of propagated samples. If the position vectors for the two satellites are correlated, then the all-on-all comparison violates the requirement that ni ¼ n0i . Hence, all-on-all solutions are not valid for correlated satellites. The method of identifying a collision between two realizations in Eqs. (22)–(25) assumes a spherical keep-out radius for the spacecraft. Such an assumption may be removed by altering these equations and improving the fidelity of the positive event criteria (Eq. (25)). For example, when given cuboid models of the two satellites, the hyper-plane separation theorem may be employed to identify a collision. de Vries and Phillion (2010) consider threedimensional models of two colliding satellites in a Monte Carlo framework, and such procedures may be used with the polynomial chaos-based methods. Alternatively, more general models for collision detection (e.g., see Lin et al., 1998) may be employed to identify a positive case. The additional computation time required is then reduced by using Eq. (22) and a sufficiently large spherical value R to filter out cases with sufficiently large separations. 3.2. PC estimation of P c and P 0s Estimating the probability of collision via PC differs slightly from the above procedure in that using a polynomial surrogate requires significantly fewer evaluations of the ODE solver. To compute P c or P 0s : 1. Use the methods described in Sections 2.1–2.1.3 with the sample generation defined by Eq. (21) to produce the ^ l ðtk ; nÞ for tk ; k ¼ 1; . . . ; K. PCE’s X ^ l ðtk ; ni Þ using Eq. 19 and 2. Generate M MC realizations X MC the set of random inputs fni gM i¼1 . 3. Compute the separation distances sðtk ; ni ; n0i Þ (Eq. (22)) for the two satellites under consideration using the realizations from step 2. 4. If needed, compute P c using Eq. (23). 5. To compute P 0s , evaluate Eqs. (24) and (25). With this procedure, the more computationally efficient polynomial evaluations replace the orbit propagators required for Monte Carlo analysis. In general, orbit propagation may take approximately 1–2 s of computer processing time per evaluation, but the polynomial evaluations (step 2) are on the scale of thousands of evaluations per second. Hence, the PC-based methods reduce the computation time needed to generate the realizations used in estimating P c and P 0s . Of course, this procedure is affected by the accuracy of the PCE. Methods exist to account for errors in the PCE

when calculating a failure probability, e.g., satellite collision, by combining PCE evaluations with additional Monte Carlo orbit propagations (Li and Xiu, 2010; Li et al., 2011). Such techniques are not required for the test cases in Section 4. Additionally, this work assumes the propagation tool is accurate and does not include any force model uncertainties. This is also an issue for Monte Carlo methods with the same random inputs to the system. Random errors in other parameters, e.g., the coefficient of drag, may be included as additional inputs to the PCE to account for the modeling errors (Jones et al., 2013). In such cases, the stochastic solution also accounts for model uncertainty at the expense of increasing d. Like the Monte Carlo methods discussed in Dolado et al. (2011), importance sampling may also be used with the PCE solution. Although M 0PC would not be reduced, such techniques will decrease the number of polynomial evaluations required to compute P c and P 0s . This is part of the technique proposed in Li et al. (2011) to mitigate for PCE approximation errors. 3.3. Summary of PC-based conjunction assessment assumptions In summary, the assumptions employed when using PC for conjunction assessment in this paper are: 1. the random inputs n are Gaussian with zero mean and identity covariance, 2. the times of interest tk ; k ¼ 1; . . . ; K, are sufficiently dense to detect a collision, 3. the black-box integrator provides sufficient accuracy to perform conjunction assessment, 4. the volume defining a positive collision is spherical, and 5. a PCE may be generated that meets the required approximation accuracy. The effects of each of these assumptions is determined by one of several design parameters, and may be mitigated through changes in the system setup. The methods for removing or reducing the effects of these assumptions were discussed in the previous sections when introduced. Additionally, although it is not commonly listed, assumption (3) above also affects the classic methods (Alfriend et al., 1999; Patera, 2001, 2003; Alfano, 2005; Chan, 2004; Alfano, 2006). 4. Numerical examples 4.1. Simulation methods This section describes the common simulation methods employed in the numerical examples presented in the following sections. Details specific to a given test case, e.g., orbit parameters, are described in later sections. The CU-TurboProp (Hill and Jones, 2009) orbit propagation suite serves as the principle simulation tool for all

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875 Table 1 Test case force models.

1867

Table 2 Test case 1 mean state (km and km/s) at t0 ¼ Feb. 3, 2016, 15:53:00 UTC.

Force

Model

Parameter

Satellite 1

Satellite 2

Gravity perturbations Sun/Moon ephemeris Geocentric terrestrial reference frame Atmospheric density Coefficient of drag (C D ) Coefficient of reflectivity (C R )

GGM02C DE421 IAU 2000AR06/2006 US Standard, 1976 2.2 1.8

X Y Z X_ Y_ Z_

35159.238 22444.088 15588.551 2.547915 0.227987 1.328323

35156.158 22441.713 15597.226 2.547674 0.228087 1.328855

4.2. Test case 1: satellites in formation tests in this paper. It provides a suite of orbit propagation capabilities using both a collection of high- and low-fidelity models of the forces acting on a satellite. Both cases consider the same forces, although with potentially different fidelities appropriate for the given orbit. A description of the force models used may be found in Table 1. As specified for each test case, the propagation tools use either a Dormand-Prince 8(7) (Prince and Dormand, 1981) or eighthorder Gauss-Jackson (Jackson, 1924; Berry and Healy, 2004) numerical integration method. The variable-step propagations using the Dormand-Prince 8(7) employ a relative tolerance of 1014. The gravity model and PCE software is implemented in Fortran, while all other computationally intensive tools are written in C. Python is then used to facilitate higher-level interaction between of the various tools. All execution time comparisons were conducted on the University of Colorado’s Janus supercomputer. Janus uses Redhat 5.5 Enterprise on 1,368 computation nodes, and each node contains two hex-core 2.8 GHz Westmere processors (12 cores per node). Software was compiled using the Intel Compiler suite (version 2011.0.013) with -O2 optimization. Tests profiling the performance of the PC-based algorithm include a comparison with the classic method presented in Alfano (2005), which is dubbed the Alfano method in this paper. The derivation of this method uses the simplifying assumptions discussed in Section 1, but with the advantage of reducing the computation cost of collision probability estimation. Each test case presented in the following sections violates two of the simplifying assumptions, thereby implying that the Alfano method may yield erroneous results. Including this method on the comparisons provides a differentiation in performance when compared to this classic technique, and is not intended to imply a shortcoming of the Alfano method other than the effect of violating the clearly defined simplifying assumptions. The implementation of the Alfano method used here is based on the presentation in Alfano (2005) and was verified by comparing to the test cases in Alfano (2009). These comparisons yielded accuracies in the current implementation that exceeded the required precision of this study by at least two orders of magnitude. In all test cases, propagation of the Gaussian uncertainty required for the classic method used the unscented transformation as presented by Van der Merwe et al. (2001).

The first test case considers two satellites in the Magnetospheric Multiscale (MMS) mission formation (Sharma and Curtis, 2005), and, in the method of solution presented, violates the assumptions of uncorrelated state vectors and nonlinear relative motion. The MMS formation includes four spacecraft in a highly elliptical orbit (eccentricities greater than 0.8) with possible collisions at the ascending and descending nodes when the orbit planes cross (Carpenter et al., 2011). Table 2 provides the initial state vectors for two simulated MMS satellites, and the position and velocity 1r uncertainties are approximately 1 m and 104 m/s, respectively, at t0 . These values are commensurate with the expected navigation uncertainties for the MMS mission at a similar point in the orbit6. Table 3 provides an overview of the conjunction at the time of closest approach (tCA ) for the mean orbit solutions. The correlations between the satellite states are described below, but violate the assumed independence of the position solutions for the two spacecraft. This test case considers two state correlation scenarios at tCA . Case 1A includes cross-satellite correlations approximately equal to 0.2, and Case 1B has correlations of 0.99. The mean solutions at t0 are the same for each case. Fig. 1 provides a histogram of the cross-satellite correlation coefficients at t0 , not tf , for both cases. The correlation coefficients for Case 1A at t0 have a maximum absolute value of 0.42, while the maximum absolute value for Case 1B is 0.93. For a true conjunction between two satellites in a constellation, as opposed to this synthetic scenario, the correlations would be a function of the observation types and geometry during the period before t0 . Cases 1A and 1B are designed to consider two possible scenarios, one a worst case, at tCA . Both test cases use the force models comparable to those prescribed for long-term studies of the MMS mission (Hughes, 2008). Given the high eccentricity and small interactions with the Earth’s aspherical gravity field, simulation of these satellites uses only a 10  10 gravity field. The area-to-mass ratio used in both drag and solar radiation pressure calculations is  1:8  103 m2/kg. For the initial propagation to tCA  300 s, the Dormand Prince 8(7) integrator is used. The Gauss-Jackson algorithm provides the

6

Personal correspondence with Dan Mattern, a.i. solutions, May 2012

1868

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

Table 3 Test case 1 conjunction information. Conjunction property

Value

tCA Distance Speed

t0 + 407315.48 s (4.71 d) 65.35 m 0.90 m/s

Fig. 1. Cross-satellite correlation coefficients at t0 for Cases 1A and 1B.

solutions for the time period of tCA  300 s with states provided every 1 s. Fig. 2 illustrates the relative dynamics of the mean solutions for the two spacecraft and P c during the time period ½tCA  300; tCA þ 300 s. The relative position is projected onto the conjunction plane with basis vectors ^ ¼ r1 ðtCA Þ  r2 ðtCA Þ ; X kr1 ðtCA Þ  r2 ðtCA Þk2 ^ ¼ ðr1 ðtCA Þ  r2 ðtCA ÞÞ  ðr_ 1 ðtCA Þ  r_ 2 ðtCA ÞÞ ; Y kðr1 ðtCA Þ  r2 ðtCA ÞÞ  ðr_ 1 ðtCA Þ  r_ 2 ðtCA ÞÞk2

ð27Þ ð28Þ

i.e., the plane perpendicular to the relative velocity at tCA , and the X component at tCA is the minimum separation of the mean solutions. The ‘linear’ solution presented in Fig. 2 is the relative trajectory based on the assumption

of linear motion. As seen in this figure, the relative dynamics of the two spacecraft in the conjunction plane is nonlinear and the relative velocity at tCA is nearly zero. Due to the latter property, the predicted relative position based on the linear assumption changes little over time. The provided P c in Fig. 2 assumes R equals 100 m. This value is relatively large, but the hard-body radius used for MMS is 120 m (Carpenter et al., 2011). To accurately assess the collision risk with the extended time period when P c P 102 , a multi-epoch analysis is required for this scenario. Fig. 3 describes the pdf of the absolute states for satellites 1 and 2 at tCA for both cases 1A and 1B. As seen in the figure, the difference in initial correlation coefficients only slightly changes the mean and standard deviation for satellite 2. The pdf for both satellites is effectively Gaussian, implying a low-degree PCE is required to approximate the stochastic solution. For a Gaussian distribution of X, a p ¼ 1 PCE exactly approximates the solution. The pdf remains Gaussian for this case due to the relatively small initial state uncertainty at t0 . PCE representations of the satellite states at the times of interest were generated with d ¼ 12; p ¼ 1; P ¼ 13, and a convergence criteria of 0.1 m. The PCE then requires M 0PC ¼ 24, requiring 48 total propagations (24 per satellite). The RSS position error of the realizations generated for the Case 1A PCE’s may be found in Fig. 4. These histograms illustrate the magnitude of the approximation error Dri (from Eq. (15)) when comparing M MC ¼ 106 Monte Carlo samples to realizations generated using the PCE with the same random inputs. The RMS errors for satellite 1 and 2 for this case were 9.4 and 12.8 cm, respectively, demonstrating that the PCE generation algorithm of Section 2.1.2 yields an approximation similar to the user-specified convergence tolerance. These approximation errors are at least an order of magnitude less than the meters scale of R, and, for this reason, the remainder of this section does not correct the computation of P 0s for the PCE approximation errors. As mentioned previously, methods exist to account for PCE errors in the computation of failure probabilities (Li and Xiu, 2010; Li et al., 2011). The resulting values of P 0s for Case 1A may be found in Fig. 5. These results use 108 Monte Carlo trials and PCE evaluations, but the random inputs used to generate the Monte Carlo trails are not the same as those used in the evaluation of the PCE’s. The abscissa is truncated at R ¼ 100 m since there is litte change in the solution performance at higher values, and, when R P 12:7 m, results indicate that P 0s P 105 . The relative precision (bottom image) quantifies the agreement between the PC- and Monte Carlo-based solutions. The estimated convergence error in a Monte Carlo-determined value of a binomial probability is (Vose, 2000, p. 149)

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ð1  P Þ Fig. 2. Relative dynamics of the Case 1 satellites in the conjunction plane rP ¼ M MC (top) and P (bottom). c

ð29Þ

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

1869

Fig. 3. Absolute position pdf and r at tCA for Cases 1A and 1B. The top row describes the position pdf for satellite 1, and the bottom row describes satellite 2.

Fig. 4. Precision of PCE-based realizations for Case 1A with the propagated solutions.

and de Vries and Phillion (2010) demonstrate that Eq. (29) yields an accurate assessment of convergence error for satellite collision probability estimation. The dashed line in the bottom plots of Figs. 5 and 6 indicates the ratio of the estimated 3rP error with the Monte Carlo determined P 0s . The PC-based solutions for Case 1A agree with the Monte Carlo results to the precision determined via the estimated convergence error. As stated in Section 1, probabilities in the range ½105 ; 1:0 are of particular interest for conjunction assessment, and M MC ¼ 108 yields a solution with an estimated 3rP accuracy of one digit or greater for such cases. As seen in the bottom plot of Fig. 5, the relative precision between the PC solution and the Monte

Fig. 5. P 0s solutions (top) and precision (bottom) for Case 1A using the methods of Alfano, PC, and Monte Carlo. The vertical lines at 12.7 m correspond to the Monte Carlo-based solution where P 0s ðRÞ  105 , and the horizontal dashed line in the lower plot indicates one digit of precision (101 ).

Carlo reference is less than 101 for a hard-body radius larger than 12.7 m. Hence, the PC-based solution for this test case provides at least one digit of accuracy when P 0s P 105 . Additionally, for R < 20 m, the Alfano method indicates a larger risk of collision, thereby increasing the chance of an unnecessary avoidance maneuver. The Case 1B results in Fig. 6 demonstrate a larger error in P 0s when using the Alfano method with highly correlated uncertainties. In this case, the risk of collision when R < 64 m is less than all decision thresholds, yet the Alfano method predicts a larger risk. In this case, there is a

1870

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

where tPC and tMC are the results for the PC- and Monte Carlo-based solutions, respectively. The computation time of each component is provided in Table 4 and, as previously stated, M 0PC is 24. These execution times include the processing over the complete 601 s interval, with average computation times computed during the generation of the data given in Fig. 5 with M MC ¼ 108 . The configuration of the PCE will effect the runtime, i.e, the value of teval changes with d and p. To generate a P 0s with an estimated 1r accuracy of one digit, Monte Carlo methods require 600 h of processing time. The PC-based methods require 13 min to perform the same comparison. 4.3. Test case 2: satellites with non-Gaussian solutions Fig. 6. P 0s solutions (top) and precision (bottom) for Case 1B using the methods of Alfano, PC, and Monte Carlo.

lower precision between the PC and Monte Carlo methods at the point where P 0s has a rapid change from 1.0 to less than 106 . The error in the PC solution likely results from the approximation error in the PCE, which has a stronger influence when P 0s is highly sensitive to sðtk ; n; n0 Þ. However, the PC methods are able to predict a near-zero P 0s for small R. A comparison of the serial computation time required for PC- and Monte Carlo-based methods of conjunction assessment for this MMS-like case may be found in Fig. 7. The figure illustrates a serial implementation, and does not account for any parallelization employed in this study. To generate the results in the figure, P 0s is assumed to be 105 , and M MC ¼ 103 ; 104 ; . . . ; 109 . The values on the abscissa are based on the assumed P 0s , the input number of samples, and Eq. (29). The ordinate provides the estimated time required to compute a solution, generate and evaluate the PCE’s (in the case of PC methods), and perform the M MC comparisons. The total computation times are estimated via   tPC ¼ tgen þ 2M 0PC tprop þ M MC teval þ tcomp   tMC ¼ M MC 2tprop þ tcomp

ð30Þ ð31Þ

The second case considers two independent objects in low-Earth orbit (LEO) with a possible collision after 48 h. For this case, the position state uncertainty during the conjunction is non-Gaussian and changes over time, thereby violating two assumptions in the Alfano method. The scenario is based on test case 7 of Alfano (2009), but with some alterations for the current paper. The propagation methods in Alfano (2009) only used two-body dynamics, so adjustments (described below) are required in the initial state to account for higher-fidelity modeling. Additionally, the initial covariance matrix is scaled by 1202 to yield skewness in the Cartesian position pdf at tCA . This yields an RSS 1r position and velocity uncertainty of 49.5 m and 0.02 m/s, respectively, for each satellite. We note that demonstrations of uncertainty propagation in the literature exhibit a more non-Gaussian solution (e.g., Jones et al., 2013), but yield a larger solution volume. In such cases, the probability of collision is very small (<105 ), which implies a low-risk scenario. Hence, such a case is not considered in this paper. Case 2 uses the higher-fidelity force models required for LEO satellites. Given the lower orbit altitude, the simulation uses a 200  200 gravity field, but represented using the cubed-sphere model to reduce computation time (Jones et al., 2010). The area-to-mass ratio used in both the drag and solar radiation pressure model is 0.1 m2/kg. This scenario only uses the Gauss–Jackson propagator, with a 10 s step size until tCA  100 s, and then a 1 s step size through tCA þ 100 s.  0 Þ, the Xðt  CA Þ vector for each satellite To generate Xðt given in Alfano (2009) is propagated back in time for 48 h to t0 using the higher-fidelity force models. This state  0 Þ used in this paper, with the then serves as the new Xðt Table 4 Test case 1 average computation times.

Fig. 7. Computation cost of serial Monte Carlo- and PC-based solutions for P c .

Process

Time (s)

Propagation of single satellite (tprop ) State comparison (tcomp ) PCE generation (tgen ) PCE evaluation (teval )

0.11 3:4  105 0.04 4:2  105

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875 Table 5 Test case 2 mean state (km and km/s) at t0 ¼ Mar. 1, 2012, 00:00:00 UTC. Parameter

Satellite 1

Satellite 2

X Y Z X_ Y_ Z_

6777.828 1085.564 210.326 0.688773 5.351702 5.387790

6780.042 1068.401 227.858 0.659875 5.357861 5.385113

Table 6 Test case 2 conjunction information. Conjunction property

Value

tCA Distance Speed

172,800 s (48 h) 3.4 m 3.20 m/s

results for both satellites provided in Table 5. Properties of the conjunction are also provided in Table 6, and indicate a closer approach between the two satellites in comparison with test case 1. The relative position (in the conjunction plane) and P c is provided in Fig. 8. The distance of the mean solutions at tCA is approximately 3.4 m with linear motion over the period tCA  50 s. Hence, this case does not violate the assumption of linearity. However, like the previous case, a multi-epoch analysis is still required. Fig. 9 illustrates the X–Y two-dimensional position pdf at tCA  50 s and tCA . The Monte Carlo results in this figure are based on 105 realizations, and exhibit skewness in the solution. A logarithmic scale is used in the color scale to more easily visualize the tails of the distribution. As seen in the figure, the 3r probability ellipsoids overlap at tCA  50 s, but the lack of intersection in the Monte Carlo-determined pdfs imply a smaller P c . As indicated by Fig. 8, P c < 105 at tCA  50 s. Additionally, the orientation of the ellipsoids change in the time period ½tCA  50; tCA . Since the solution is both non-Gaussian

Fig. 8. Relative dynamics of the Case 2 satellites in the conjunction plane (top) and P c (bottom).

1871

and the uncertainty changes during the conjunction, two of the common collision probability estimation assumptions are violated in this case. Finally, the right-most plots compare the Monte Carlo realizations to the propagated Gaussian pdf. The Gaussian approximation predicts a more diffuse solution at tCA , implying a lower instantaneous probability of collision. Table 7 provides the P c at these two times based on a 107 sample Monte Carlo analysis and the same number of samples generated from the propagated Gaussian pdf. The table confirms the qualitative assessment of Fig. 9. PCE representations of the satellite states at the times of interest are generated with d ¼ 6; p ¼ 3; P ¼ 84, and a PCE generation convergence tolerance of 0.1 m. The PCE for satellite 1 requires M 0PC;1 ¼ 116, while satellite 2 requires M 0PC;2 = 148 training samples. Note the extra subscript denotes satellite one or two. Although it is not readily apparent in Fig. 9, the Y and Z Cartesian component standard deviations for satellite 2 are slightly larger than those for satellite 1. To approximate the stochastic solution with the increased solution distribution, the satellite 2 PCE solution requires additional samples for its generation. The RSS position error of the realizations generated for the Case 2 PCE’s are given in Fig. 10. These histograms summarize the comparison of 106 Monte Carlo samples to realizations computed using the PCE with the same random inputs, and illustrate that the PCE generation algorithm yields an estimate with accuracies approximately equal to the generation algorithm convergence tolerance. Fig. 11 provides the estimated P 0s for Case 2 using both the one-to-one and all-on-all sampling methods. The reference for all solutions is the one-to-one Monte Carlo analysis using M MC ¼ 108 trials. The provided one-to-one PC estimate also used 108 PCE evaluations. The all-on-all solutions use 10,000 trials per satellite, but reuse samples in Eq. (22) to reduce the number of propagations. This yields 108 comparisons in the all-to-all analysis. Both the PC oneto-one and all-to-all solutions agree with the reference to the estimated 3r convergence accuracy. Additionally, the PC solution provides at least one digit of accuracy when P 0s P 105 . Since two fundamental assumptions used in the derivation of the Alfano method are not satisfied in this case, it fails to provide an accurate collision probability in the critical decision regions of P 0s 2 ½105 ; 103 . Hence, the Alfano method can cause a failed collision detection for some values of R. Fig. 12 presents a comparison of the expected computation time required for Case 2. Average execution times for each of the principle operations is provided in Table 8. The method for generating the top result is the same as that of Fig. 7, but the bottom image of Fig. 12 adjusts the number of samples for all-on-all methods. For this latter plot, the estimated execution times are

pffiffiffiffiffiffiffiffiffiffiffi t0PC ¼ tgen;1 þ tgen;1 þ M 0PC;1 þ M 0PC;2 tprop þ 2teval M MC þ tcomp M MC pffiffiffiffiffiffiffiffiffiffiffi t0MC ¼ 2tprop M MC þ tcomp M MC ;

ð32Þ ð33Þ

1872

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

Fig. 9. Case 2 satellite pdf in the X–Y plane at tCA  50 s (top) and tCA (bottom). The leftmost column provides the 3r position ellipsoids, the center is the Monte Carlo-determined pdf (logarithmic color scale), and the right column is the combination of the previous two.

Table 7 Monte Carlo- and Gaussian-based P c ðR ¼ 100 mÞ. Time tCA  50 s tCA

Gaussian 5

8.710 2.6103

Monte Carlo 4.8106 3.0103

Fig. 10. Precision of PCE-based realizations for Case 2 with the propagated solutions.

pffiffiffiffiffiffiffiffiffi where the M PC terms result from performing the all-onall analysis instead of a one-to-one. Like Case 1, the presented total computation times assume no parallelization

Fig. 11. P 0s solutions for Case 2 using the methods of Alfano, PC, and Monte Carlo (top), relative accuracy of the one-to-one PC solution with the comparable Monte Carlo solution (middle), and the relative precision of the all-on-all PC and Monte Carlo solutions when compared to the oneto-one Monte Carlo solution. The vertical lines at 4 m correspond to the Monte Carlo-based solution where P 0s ðRÞ  105 , and the horizontal dashed line in the lower plots indicate one digit of precision (101 ).

and incorporate operations over the full conjunction period (401 s). As seen in the figure, Monte Carlo methods demand almost 6,000 h to generate a P 0s that is accurate to

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

Fig. 12. Computation cost of serial Monte Carlo- and PC-based solutions for P 0s when using one-to-one (top) and all-on-all (bottom) sampling.

Table 8 Test case 2 average computation times. Process

Time (s)

Propagation of single satellite (tprop ) State comparison (tcomp ) PCE generation (tgen;1 =tgen;2 ) PCE evaluation (teval )

1.01 2:5  105 0.013/0.016 3:8  105

approximately one digit with a one-to-one analysis, but only 112.8 min are needed with an all-on-all sampling. The PC-based methods require approximately 22 and 8.4 min for the one-to-one and all-on-all methods, respectively. Since the time required for a one-to-one analysis using PC is only three times that of an all-to-all comparison, the latter method is not required to make the problem tractable. Initial inspection of the bottom plot in Fig. 12 implies the execution time of the PC-based methods approach that of the Monte Carlo methods, but this is only an artifact of the logarithmic scale. Upon deriving t0MC  t0PC from Eqs. (32) and (33), the computation time increases faster for Monte Carlo methods since tprop > teval .

4.4. Discussion of execution time performance As indicated in the previous sections, the single-processor implementations of the PC-based methods require minutes to perform the required analysis when P 0s  105 . This is dominated by the time required to perform the M MC comparisons between the realizations. Although it provides improved solution efficiency for cases that require a particle-like formulation, such a time is intractable when considering tens of thousands of comparisons (such as in the SSA problem). Several methods exist to reduce this computation time, including: (i) parallelization of the PCE evaluation and state comparisons, (ii) importance sampling to reduce

1873

M MC , and (iii) a Galerkin projection of the absolute-state PCEs to form a stochastic solution that provides sðtk ; n; n0 Þ. Both the Monte Carlo- and PC-based methods may leverage parallel implementations and importance sampling, which will reduce the time required to compute P 0s . As described in Section 2.1.3, parallelization of the PC methods may leverage already existing software to reduce software development time. The simple, ‘embarrassingly parallel’ mathematical formulation of Eq. (19) also allows for implementation on a GPU, which is a low-cost option when compared to multi-core desktop and supercomputing systems. However, such an implementation must be tested to quantify the reduction in wall-clock time. When given the polynomial surrogate, the same importance sampling methods used with Monte Carlo (e.g., Dolado et al., 2011) may be employed with polynomial evaluations substituted for the orbit propagations. This also reduces the number of comparisons required and the total computation time. In the case of no correlation between objects, the two absolute-state PCEs may be used with a Galerkin projection to form a PC-based approximation of the relative separation sðtk ; n; n0 Þ. The accuracy and computation cost of such a solution has not yet been demonstrated. Since s (Eq. (22)) is a nonlinear equation of the two absolute states r1 and r2 , the resulting PCE will require a larger p and numerical approximation methods for the generation of ca from the original PCEs (Le Maitre and Knio, 2010). Such a process is a one-time cost at each time tk , and would increase tgen . However, each realization of the new PCE yields sðtk ; n; n0 Þ, thereby eliminating the M MC evaluations of Eq. (22) and reducing tcomp . Given the potential reduction in the variable-cost computation time, such work will be considered in future development of PC-based methods of conjunction assessment. 5. Conclusions This paper demonstrated the use of polynomial chaos (PC) for conjunction assessment. The PC-based methods provide a polynomial surrogate, called a polynomial chaos expansion (PCE), to map random inputs at the initial time to generate a realized state at a given time in the future. This allows for Monte-Carlo-like estimation of collision probabilities using polynomial evaluations instead of computationally expensive orbit propagations. The discussion provided the details required to implement the PC methods, which included a multi-epoch analysis of the conjunction to assess the survival probability over the full period of conjunction. This eliminates the need to accurately solve for the time of closest approach when a conjunction start and end time are known. Numerical demonstrations of the PC methods considered two test cases designed to demonstrate the use of the polynomial surrogates in situations where classic methods of estimating the probability of collision are not applicable. This includes a test with nonlinear relative motion and

1874

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875

coupled state uncertainties, and another case with a nonGaussian a posteriori probability density function that changed over the period of conjunction. For both of these test cases, the PCE’s provided a probability of collision that agrees with the Monte Carlo baseline solutions to the convergence accuracy of the reference. Additionally, each PC-based analysis required minutes to generate an approximation of the probability of collision, whereas a serial implementation of Monte Carlo would take tens of minutes to hours to yield the same result. As mentioned in the introduction, this work assumed an a priori identification of the conjunction and the period of possible collision. Follow-up work includes the development of PC-based techniques that address other problems in conjunction assessment, e.g., filtering of a space object catalog to identify a potential collision for further study. The PC formulation also allows for the inclusion of stochastic force model parameters, which will be introduced in future research. Finally, the polynomial surrogates enables post-maneuver collision assessment, especially when the a priori description of the maneuver uncertainty is non-Gaussian, which aids in the design of a conjunction maneuver. Acknowledgments This work was partially funded by the NASA Goddard Space Flight Center and a.i. solutions. The authors wish to thank Russell Carpenter (NASA/GSFC), Geoff Wawrzyniak, and Dan Mattern (a.i. solutions) for helping to understand the MMS mission, which served as the basis for the first test case presented in the paper. Alireza Doostan was also supported by the United States Department of Energy under Advanced Scientific Computing Research Early Career Research Award DE-SC0006402. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver and the National Center for Atmospheric Research. References Alfano, S., 1994. Determining satellite close approaches, Part II. The Journal of the Astronautical Sciences 42, 143–152. Alfano, S., 2005. A numerical implementation of spherical object collision probability. The Journal of the Astronautical Sciences 53, 103–109. Alfano, S., 2006. Addressing nonlinear relative motion for spacecraft collision probability. In: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Keystone, CO. Alfano, S., 2009. Satellite conjunction monte carlo analysis. In: 19th Annual AAS/AIAA Spaceflight Mechanics Meeting, Savannah, GA, 2009. Alfano, S., 2012. Toroidal path filter for orbital conjunction screening. Celestial Mechanics and Dynamical Astronomy 113, 321–334. Alfano, S., Greer, M.L., 2003. Determining if two solid ellipsoids intersect. Journal of Guidance, Control, and Dynamics 26, 106–110.

Alfano, S., Negron Jr., D., 1993. Determining satellite close approaches. The Journal of the Astronautical Sciences 41, 217–225. Alfriend, K.T., Akella, M.R., Frisbee, J., Foster, J.L., Lee, D.J., Wilkins, M., 1999. Probability of collision error analysis. Space Debris 1, 21–35. Armellin, R., Morselli, A., Di Lizia, P., Lavagna, M., 2012. Rigorous computation of orbital conjunctions. Advances in Space Research 50, 527–538. Babusˇka, I., Tempone, R., Zouraris, G.E., 2005. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM Journal on Numerical Analysis 42, 800–825. Berry, M.M., Healy, L.M., 2004. Implementation of Gauss–Jackson integration for orbit propagation. Journal of Astronautical Sciences 52, 331–357. Carpenter, J.R., 2002. Partially decentralized control architectures for satellite formations. In: AIAA Guidance, Navigation, and Control Conference and Exhibit, Monterey, CA. Carpenter, J.R., Markley, F.L., Alfriend, K.T., Wright, C., Arcido, J., 2011. Sequential probability ratio test for collision avoidance maneuver decisions based on a bank of norm-inequality-constrained epochstate filters. In: AAS/AIAA Astrodynamics Specialist Conference, Girdwood, AK. Chan, K., 2004. Short-term vs long-term spacecraft encounters. In: AIAA/ AAS Astrodynamics Specialist Conference and Exhibit, Providence, RI. de Vries, W.H., Phillion, D.W., 2010. Monte Carlo method for collision probability calculations using 3D satellite models. In: Proceedings of the Advanced Maui Optical and Space Surveillance Technologies Conference, Wailea, Maui, Hawaii. DeMars, K.J., Jah, M.K., Giza, D., Kelecy, T., 2009. Orbit determination performance improvements for high area-to-mass ratio space object tracking using an adaptive Gaussian mixtures estimation algorithm. In: 21st International Symposium on Space Flight Dynamics, Toulose, France. Dolado, J.C., Legendre, P., Garmier, R., Revelin, B., Pena, X., 2011. Satellite collision probability computation for long term encounters. In: AAS/AIAA Astrodynamics Specialist Conference, Girdwood, AK. Doostan, A., Ghanem, R.G., Red-Horse, J., 2007. Stochastic model reduction for chaos representations. Computational Methods in Applied Mechanical Engineering 196, 3951–3966. Doostan, A., Iaccarino, G., 2009. A least-squares approximation of partial differential equations with high-dimensional random inputs. Journal of Computational Physics 228, 4332–4345. Doostan, A., Owhadi, H., 2011. A non-adapted sparse approximation of PDEs with stochastic inputs. Journal of Computational Physics 230, 3015–3034. Doostan, A., Validi, A., Iaccarino, G., 2013. Non-intrusive low-rank separated approximation of high-dimensional stochastic models. Computation Methods in Applied Mechanical Engineering 263, 42–55. Fujimoto, K., Scheeres, D.J., 2012. Non-linear propagation of uncertainty with non-conservative effects. In: 22nd Annual AAS/AIAA Spaceflight Mechanics Meeting, Charleston, SC. Fujimoto, K., Scheeres, D.J., Alfriend, K.T., 2012. Analytical nonlinear propagation of uncertainty in the two-body problem. Journal of Guidance, Control, and Dynamics 35, 497–509. Garmier, R., Dolado, J., Pena, X., Revelin, B., Legendre, P., 2011. Collision risk assessment for multiple encounters. In: European Space Surveillance Conference, Madrid, Spain. Ghanem, R.G., 1999. Ingredients for a general purpose stochastic finite elements implementation. Computer Methods in Applied Mechanics and Engineering 168, 19–34. Ghanem, R.G., Dham, S., 1998. Stochastic finite element analysis for multiphase flow in heterogeneous porous media. Transport in Porous Media 32, 239–262. Ghanem, R.G., Red-Horse, J., 1999. Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach. Physica D: Nonlinear Phenomena 133, 137–144. Ghanem, R.G., Spanos, P.D., 1991. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag, New York.

B.A. Jones, A. Doostan / Advances in Space Research 52 (2013) 1860–1875 Ghrist, R.W., Plakalovic, D., 2012. Impact of non-Gaussian error volumes on conjunction assessment risk analysis. In: AIAA/AAS Astrodynamics Specialist Conference, Minneapolis, MN. Hill, K., Jones, B.A., 2009. TurboProp Version 4.0. Colorado Center for Astrodynamics Research. Hoots, F.R., Crawford, L.L., Roehrich, R.L., 1984. An analytic method to determine future close approaches between satellites. Celestial Mechanics 33, 143–158. Horwood, J.T., Aragon, N.D., Poore, A.B., 2011. Gaussian sum filters for space surveillance: theory and simulations. Journal of Guidance, Control, and Dynamics 34, 1839–1851. Hosder, S., Walters, R.W., Perez, R., 2006. A non-intrusive polynomial chaos method for uncertainty propagation in CFD simulations. In: 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV. Hughes, S.P., 2008. Formation design and sensitivity analysis for the magnetospheric multiscale mission (MMS). In: AIAA/AAS Astrodynamics Specialist Conference, Honolulu, Hawaii. Jackson, J., 1924. Note on the numerical integration of d 2 x=dt2 ¼ f ðx; tÞ. Monthly Notices of the Royal Astronomical Society 84, 602–606. Jia, B., Cai, S., Cheng, Y., Xin, M., 2012. Stochastic collocation method for uncertainty propagation. In: AIAA/AAS Astrodynamics Specialist Conference, Minneapolis, MN. Jones, B.A., Born, G.H., Beylkin, G., 2010. Comparisons of the cubedsphere gravity model with the spherical harmonics. Journal of Guidance, Control, and Dynamics 33, 415–425. Jones, B.A., Doostan, A., Born, G.H., 2012. Conjunction assessment using polynomial chaos expansions. In: International Symposium on Space Flight Dynamics, Pasadena, CA. Jones, B.A., Doostan, A., Born, G.H., 2013. Nonlinear propagation of orbit uncertainty using non-intrusive polynomial chaos. Journal of Guidance, Control, and Dynamics 36, 430–444. Le Maitre, O.P., Knio, O., 2010. Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics. Springer. Le Maitre, O.P., Najm, H.N., Ghanem, R.G., Knio, O.M., 2004. Multiresolution analysis of Wiener-type uncertainty propagation schemes. Journal of Computational Physics 197, 502–531. Li, J., Li, J., Xiu, D., 2011. A efficient surrogate-based method for computing rare failure probability. Journal of Computational Physics 230, 8683–8697. Li, J., Xiu, D., 2010. Evaluation of failure probability via surrogate models. Journal of Computational Physics 229, 8966–8980. Lin, M.C., Gottschalk, S., 1998. Collision detection between geometric models: a survey. In: Proceedings of the IMA Conference on Mathematics of Surfaces, pp. 602–608. Mercurio, M., Singla, P., Patra, A., 2012. A hierarchical tree code based approach for efficient conjunction analysis. In: AIAA/AAS Astrodynamics Specialist Conference, Minneapolis, MN. National Research Council, 2012. Continuing Kepler’s Quest: Assessing Air Force Space Command’s Astrodynamics Standards, The National Academies Press. Park, R.S., Scheeres, D.J., 1367. Nonlinear mapping of Gaussian statistics: theory and applications to spacecraft trajectory design. Journal of Guidance, Control, and Dynamics 29, 1367–1375.

1875

Patera, R.P., 2001. General method for calculating satellite collision probability. Journal of Guidance, Control, and Dynamics 24, 716– 722. Patera, R.P., 2003. Satellite collision probability for nonlinear relative motion. Journal of Guidance, Control, and Dynamics 26, 728–733. Prince, R.J., Dormand, J.R., 1981. High order embedded Runge–Kutta formulae. Journal of Computational and Applied Mathematics 7, 67– 75. Sabol, C., Binz, C., Segerman, A., Roe, K., Schumacher, Jr., P.W., 2011. Probability of collision with special perturbations dynamics using the monte carlo method. In: AAS/AIAA Astrodynamics Specialist Conference, Girdwood, AK. Sabol, C., Sukut, T., Hill, K., Alfriend, K.T., Wright, B., Li, Y., Schumacher, Jr., P.W., 2010. Linearized orbit covariance generation and propagation analysis via simple monte carlo simulations. In: 20th Annual AAS/AIAA Spaceflight Mechanics Meeting, San Diego, CA. Sharma, A.S., Curtis, S.A., 2005. Magnetospheric multiscale mission. In: Nonequilibrium Phenomena in Plasmas. Astrophysics and Space Science, Library, vol. 321. Springer, The Netherlands, pp. 179–195. Soize, C., Ghanem, R., 2004. Physical systems with random uncertainties: chaos representations with arbitrary probability measure. SIAM Journal of Scientific Computing 26, 395–410. Valli, M., Armellin, R., Di Lizia, P., Lavagna, M.R., 2013. Nonlinear mapping of uncertainties in celestial mechanics. Journal of Guidance, Control, and Dynamics 36, 48–63. Van der Merwe, R., Wan, E.A., 2001. The square-root unscented Kalman filter for state and parameter-estimation. In: 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing, Salt Lake City, Utah, pp. 3461–3464. Vose, D., 2000. Risk Analysis: A Quantitative Guide, Second ed. John Wiley and Sons, Ltd., Chichester, England. Wan, X., Karniadakis, G.E., 2005. An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. Journal of Computational Physics 209, 617–642. Wiener, N., 1938. The homogeneous chaos. American Journal of Mathematics 60, 897–936. Woodburn, J., Coppola, V., Stoner, F., 2009. A description of filters for minimizing the time required for orbital conjunction computations. In: AAS/AIAA Astrodynamics Specialist Conference, Pittsburgh, PA. Woodburn, J., Dichmann, D., 1997. Determination of close approaches for constellations of satellites. In: IAF International Workshop on Mission Design and Implementation of Satellite Constellations, Toulouse, France. Xiu, D., 2009. Fast numerical methods for stochastic computations: a review. Communications in Computational Physics 5, 242–272. Xiu, D., 2010. Numerical Methods for Stochastic Computations: A Spectral Method Approach. Princeton University Press, Princeton, NJ. Xiu, D., Karniadakis, G.E., 2002. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM Journal of Scientific Computing 24, 619–644. Xiu, D., Karniadakis, G.E., 2003. Modeling uncertainty in flow simulations via generalized polynomial chaos. Journal of Computational Physics 187, 137–167.