Tarek A. Elgohary University of Central Florida Internet of Things (IoT) Summit at RWW2019

1. A New Method for Computing Orbital Probability of Collision and Autonomous Space-based Orbit Estimation with Multiple Agent Nodes. Tarek A. Elgohary University of Central Florida Internet of Things (IoT) Summit at RWW2019 January 21st, 2019.

2. Motivation https://orbitaldebris.jsc.nasa.gov/

4. Inertial EOM for m1: Inertial EOM for m2: Since Then the relative EOM for m2 relative to m1 is N O The Two-Body Problem • Two-Body EOM

5. SSA & UQ • SSA (ESA): Space Weather, Near Earth Objects, and Space Surveillance and Tracking • NASA: NEOWISE, Orbital Debris at JSC, Active Debris Removal. • US-AF: Catalog Maintenance, Correlations and Collisions. • Uncertainty characterization is a major challenge for SSA • Catalog Maintenance • Uncorrelated Tracks • Space Object Characterization • Collisions • Standard methods for uncertainty quantification consider perturbations to be random and uncorrelated

6. Existing Methods: UQ Numerically solving the Fokker-Planck equation Monte Carlo Simulations Gaussian Mixture Models Polynomial Chaos Expansions

7. Existing Methods: CA • STM based propagation of uncertainty • Combined PDF methods • Foster’s Method • Chan’s Method • Patera’s Method • Alfano’s Method • All these methods rely on the assumption that the distribution has remained Gaussian • Essentially a very short time of conjunction, rectilinear motion. • May no apply for longer conjunctions, e.g. GEO.

8. Parameter Uncertainty and Liouville's theorem The Fokker–Planck equation is normally used to characterize the time evolution of Probability Density If the standard Wiener process (stochastic noise) can be replaced with parameterized uncertainty then Liouville's theorem can be substituted This allows for the correlation of initial probability with that at any time along the trajectory. Fokker–Planck Equation: Liouville's Theorem: The classic Method of Characteristics Solution Process is one possibility to solve Liouville’s Eq. However we propose a novel Inverse Method of Characteristics

9. Liouville Theorem Application Utilizing Liouville's theorem allows for a-priori probability information to improve posterior uncertainty estimates Example: 1,000 Sample Monte Carlo and kernel density estimation for a duffing oscillator The Liouville version provides significantly better results with the same number of function evaluations

10. Orthogonal Probability Approximation Problem: Known a-priori PDF, with an unknown non-Gaussian posterior PDF Concept: Use orthogonal approximation in combination with Liouville's theorem to build high accuracy representations of posterior PDFs

11. Orthogonal Probability Approximation Method: Use forward propagation to bound region of interest

12. Orthogonal Probability Approximation Method: Use forward propagation to bound region of interest Populate ROI with evaluation nodes

13. Orthogonal Probability Approximation Method: Use forward propagation to bound region of interest Populate ROI with evaluation nodes Back propagate and use Liouville to obtain posterior probability value

14. Orthogonal Probability Approximation Method: Use forward propagation to bound region of interest Populate ROI with evaluation nodes Back propagate and use Liouville to obtain posterior probability value Generate approximation for posterior PDF

15. PDF Approximation and Integration Currently ND Chebyshev polynomials are being used for approximations with standard Clenshaw and Curtis integration scheme Hermite polynomials with spares gridding are being investigated for higher dimensions to mitigate the curse of dimensionality

16. Linear Validation A simple linear oscillator with an initial Gaussian distribution was propagated using analytic linear error theory for validation of OPA

17. Linear Validation The result of Gaussian linear error propagation theory. With order 40 the approximate PDF’s accuracy matched to 1e-6, similar to the specified integration tolerance (1e-6).

18. PDF Approximation Results Total Probability P = 1.00000000 Note: Log-like scaling, 0=0 & 1=1

19. 2-DoF Undamped Duffing Oscillators Propagate extreme probability bound to TOI to identify ROI Populate 2D cosine grid Back-propagate to t0 and get probability for a-priori PDF Integrate subsequent dimensions to produce marginalized probability

20. Time Domain Approximation

21. Conjunction of Duffing Oscillators • To validate using OPA probability of collision a pair of undamped Duffing oscillators with no periodic excitation are chosen. • This provides a nonlinear two dimensional state space to test OPA.

22. Conjunction ROI Propagate extreme probability bound to TOI for both objects to identify overlapping ROI Consequential Dimensions: Dimensions where the analysis is performed. Spatial dimensions in the case of CA Marginal Dimensions: Dimensions that are marginalized into the consequential dimensions. The PDF approximation can be limited to the overlapping consequential dimensions, this reduces required computation, and if there is no intersection then the analysis is complete

23. Multiple Marginalized PDFs

24. Object Probability of Collision Once the marginalized probability is available for each object the probability of collision can be computed By switching which object is treated as the primary object can check numerical approximation convergence

25. Probability of Collision Results

26. Probability of Collision Validation Monte Carlo analysis was used to validate the performance of OPA For each possible conjunction case one million samples are propagated for the two oscillators For each sample of the first oscillator every sample of the second was checked to see if it was < away The average percentage of this result for each of the first oscillators samples is taken as the probability of collision

27. Monte Carlo Validation Results

28. Planar Orbit PDF Approximation

29. Planar Orbital Conjunction Currently we have extended OPA to a planar orbit problem Two arbitrary RSOs are set up to have a near collision Currently using only two-body acceleration so that the problem remains planar

30. Conjunction ROI In Higher Dimensions To identify the extremal probability regions in higher dimensions Monte Carlo like forward propagation is performed with initial state that have an extremely low probability as defined be the a-priori distribution Once these state have been propagated to the time of interest their overlap in the consequential dimensions determines the region of interest

31. Planar Orbital Conjunction: PDF Extending the previously described steps to four dimensions makes it possible to identify the region of possible intersection. The uncertainty in the velocity can then be marginalized generate functional approximations for the RSO’s PDFs

32. Planar Orbital Conjunction: Probability of Collision These approximations for PDFs allow for the computation of the local probability of collision and integrating that probability give the total probability We are currently running validation cases for the orbital problem with Monte Carlo

33. Orbital Conjunction: Velocity Only Approach The standard procedure for OPA is followed until the region of interest is identified and the grid of evaluation nodes is placed Each evaluation node is back propagated to the initial point using a BVP solver, such as a Lambert solver or the Method of Particular Solutions The velocity value for the initial state is compared to the a-priori PDF to obtain a probability value

34. Orbital Conjunction vs. MC Velocity only OPA was used to compute the PC for a pair of RSOs with various conjunction parameters The same analysis was computed using MC

35. OPA Updates: ROI Alignment To reduce the number of evaluation nodes that are end up evaluating negligible probabilities the output of the LMC evaluation used to find the ROI can be used to find a pseudo mean and covariance. Eigen analysis of these pseudo values can be used to realign the ROI so that it better contains the PDF.

36. OPA Updates: ROI Segmentation One issue with standard Clenshaw–Curtis quadrature is that the cosine nodes used for the approximation are least dense in the center of the approximation region. Using a pair of cosine samples can greatly improve the distribution of nodes with respect to the PDF. Generating the cosine polynomials so that they meet at the nominal result can improve the result further

37. OPA with Alternate Quadrature One major way to address the dimensionality issues of OPA is to replace Clenshaw–Curtis quadrature in the marginal dimensions with another method of quadrature This allows for the cost reductions of other quadrature methods to compute the probability contributions of the marginal dimensions, while maintaining the advantage of the Clenshaw–Curtis functional approximation in the consequential dimensions

38. OPA with Smolyak Quadrature Another alternate method that can be utilized is Smolyak sparse-grid quadrature Smolyak sparse-grid quadrature uses a sparse tensor product to produce a set of evaluation nodes

39. 3D Cosine Grid vs Smolyak Grid

40. Cosine Grid vs Smolyak Grid Node Cost The cost advantage of Smolyak sparse-grid quadrature increases substantially as the dimension number increases

41. Planar Orbit PDF: Sparse Grid Single RSO considering a planar orbit (four dimensions). The resulting PDF and the cumulative probability

42. 6-DoF Approximation Results The enhancements described were combined to allow for approximation of the PDF for 6-DoF orbits We are currently extending this approximation for POC analysis to 6+ Dimensions

43. Semi-Stochastic OPA One alternate method of quadrature is to use LMC to find the probability contribution of the marginal dimensions Each evaluation node in the consequential dimensions is a separate LMC computation This method can also account for stochastic uncertainty that can not be converted to a parametric dimension Single RSO considering a planar orbit with drag (five dimensions)

44. Non-Conservative Systems The probability density is no longer conserved when considering non-conservative systems. EG. Damped Linear Oscillator.

45. Non-Conservative Systems However, the change in the probability density is scaled with the change in the volume at each evaluation node. This can be computed by comparing the simplectic volume of the posterior grid to the deformed grid at the initial time.

46. Evaluation Node Volume Change Comparing the change in the volume around the evaluation nodes at the final time to the volume of the distorted nodal positions at the initial time The volume change is estimated by using change in distance between neighboring nodes In 2D the shoelace formula can be used. For higher dimension, the nodes can be considered to be the vertices on nD simplexes 2D nD

47. Non-Conservative Systems Preliminary results for damped Duffing Oscillator PDF approximation.

48. ARTISAN: Autonomous Real-Time Identification with Space-based Agent Nodes

50. Satellite Attitude and Angular Velocity Nonlinear Control Using Lyapunov Control Feedback Use quaternions as the form of attitude due to their non-singularity while satisfying the norm constraint that govern them. The candidate Lyapunov function for the nonlinear control is given as Where is vector component of the quaternion attitude error, whose dynamics are governed through , and is the angular velocity error. This results in a control law as follows Where is the current angular velocity, is the desired angular velocity, and and are control gains.