Intelligent Selection of NEO Deflection Strategies Under Uncertainty

This paper presents an Intelligent Decision Support System (IDSS) that can automatically assess the suitable robust deflection strategies to respond to an asteroid impact scenario. The input to the IDSS is the warning time, the orbital parameters and mass of the asteroid and the corresponding uncertainties. The output is the deflection strategies that are more likely to o ff er a successful deflection. Both aleatory and epistemic uncertainties on ephemerides and physical properties of the asteroid are considered. The training data set is produced by generating thousands of virtual impactors, sampled from the current distribution of Near Earth Objects (NEO). For each virtual impactor we perform a robust optimisation, under mixed aleatory / epistemic uncertainties, of the deflection scenario with di ff erent deflection strategies. The robust performance indices is considered by the deflection e ff ectiveness, which is quantified by Probability of Collision post deflection. The IDSS is based on a combination Dempster-Shafer theory of evidence and a Random Forest classifier that is trained on the data set of virtual impactors and deflection scenarios. Five deflection strategies are modelled and included in the IDSS: Nuclear Explosion, Kinetic Impactor, Laser Ablation, Gravity Tractor and Ion Beam Shepherd. Simulation results suggest that the proposed decision support system can quickly provide robust decisions on which deflection strategies are to be chosen to respond to a NEO impact scenario. Once trained the IDSS does not require re-running expensive simulations to make decisions on which deflection strategies are to be used and is, therefore, suitable for the rapid pre-screening or reassessment of deflection options. 23 of a deflection action has to be re-assessed multiple times given the level of uncertainty and maturity of the deflection technology and 24 the knowledge of the target NEO. Thus making optimal decisions, especially in the preliminary phase, requires a proper treatment of 25 both aleatory and epistemic uncertainty and a way to quickly assess the appropriateness of a technology under uncertainty. 26 This paper proposes an Intelligent Decision Support System (IDSS) based on a combination of ML and Dempster-Sha ff er theory 27 of evidence (DSt) to handle epistemic uncertainty. A ML model is trained on a data-set of virtual impactors, deflection scenarios 28 and deflection technologies with associated uncertainty. Once trained, the ML model works like an oracle that can be interrogated 29 multiple times to assess which deflection methods is to be used in response to a threatening scenario. In this context, Machine 30 Learning provides a way to make robust decisions under uncertainty without the need of multiple runs of expensive mission and 31 system design optimisations. The ML model encapsulates di ff erent processes that concur to assess which deflection strategies are 32 likely to o ff er a successful deflection: computation of the deflection action, transfer trajectory definition, quantification of system and 33 NEO uncertainty, optimisation of the deflection mission. DSt adds a further important layer to the decision process by quantifying the 34 uncertainty deriving from a lack of knowledge on the NEO, the deflection action and the deflection technology. The use of Machine 35 Learning has found a growing range of application in the space sector in the last decade. However, its application in the field of 36 Planetary Defense (PD) is still limited. Recent examples include the selection of asteroid deflection strategy Nesvold et al. a 37 PD Resource Discovery Engine Bambacus et al. and the estimation of momentum transfer factor β from experimentsRaskin 38 et al. Meanwhile, a decision support system for space tra ffi c management was proposed by the authors to provide operators with automatic collision avoidance capabilities The IDSS in this paper is conceptually 40 derived from


Introduction
1 Asteroid impact poses a major threaten to all life forms on Earth. Several serious impact events through history, from the Chixulub

48
This section will briefly introduce the proposed methodology and the approach to uncertainty quantification and propagation. 49 Then it will present the models used to compute both impulsive and slow push deflections and finally it will describe the approach we propose to compute robust and optimal deflection solutions. sources on the physical properties of the NEO and the technologies to be used. In this paper, as a preliminary demonstration, we 56 considered only the diameter of the NEO, with associated uncertainty, and the uncertainty on its ephemerides. The output of each 57 classifier is a label that says if the corresponding method can achieve a successful deflection or not. Thus the overall output of the 58 IDSS is a number of methods that are expected to be successfully applicable to the input scenario.
Let λ ae and λ H be two parameter vectors defining the uncertainty models ∆ae(λ ae ) and ∆H(λ H ), the expected values ae = µ T ae and H = µ T H . Then,ae andH under mixed aleatory and epistemic uncertainties are expressed as: with λ ae and λ H subject to the conditions λ ae ∈ [λ ae , λ ae ] and λ H ∈ [λ H , λ H ]. In the following the values ∆ae and ∆H are drawn from 71 two normal distributions N(0, σ T ae ) and N(0, σ T H ) respectively, thus λ ae = σ T ae and λ H = σ T H . The covariance matrix of the orbital 72 elements Σ ae at the initial state is defined as a diagonal matrix The effect of an impulsive change in the velocity of the asteroid induces a variation of its orbital elements. According to the 91 method proposed in Ref. Vasile & Colombo (2008), the deflection position on the B-plane can be calculated analytically by rewriting 92 δv in a tangential, normal, out-of-plane reference frame The variation of Keplerian elements at the time of deflection t d caused by δv can be calculated by due to the change of semi-major axis, the change in the mean anomaly at the collision time t MOID is given by therefore, the total variation in the mean anomaly δM between the unperturbed and the deflected orbit is The position of the deflected asteroid with respect to the undeflected one at the true anomaly θ MOID along the orbit of the undeflected where δr = δx r , δy θ , δz h T is the deflection vector at the collision time t MOID in a radial, transverse, out-of-plane reference 99 frame attached to the undeflected asteroid. In short, the deflection vector δr (δv, t d ; ae, ∆ae) at collision time t MOID after the impulsive 100 strategy can be expressed as where δae = δã, δẽ, δĩ, δΩ, δω, δM T = δae+∆ae indicates the total variation of orbital elements at t MOID , ∆ae indicates the uncertainty 102 of orbital elements, and A MOID indicates the transition matrix extracted from Eq.(10). In the general case of slow-push strategies, the variation of orbital elements is computed by numerical integration of Gauss order of magnitude lower than numerical integration. In the reminder of this section we will give a brief introduction of FPET. The variation of the orbital elements is obtained by integrating Gauss equations in non-singular equinoctial elements: where L is the true longitude, B = 1 +P 2 1 +P 2 2 , Φ(L) = 1 +P 1 sinL +P 2 cosL, and 114P 1 =ẽ sin(Ω +ω) P 2 =ẽ cos(Ω +ω) The components of thrust vector is formed by modulus ϵ, azimuth α and elevation β of the thrust acceleration in the radial-transverse If one assumes that the modulus of the perturbing acceleration is small compared to the local gravitational acceleration, then By substituting Eq.14 and Eq.15 into Eq.12 , the system of equations in the longitude L can be summarized in a vector form as whereẼ = [ã,P 1 ,P 2 ,Q 1 ,Q 2 ] T . The solution of Eq.16 can be expanded to the first order in the perturbing parameter ϵ whereL =L 0 + ∆L. By substituting Eq.17 into Eq.16 and expanding the right hand side in Taylor series with respect to ϵ, collecting 120 the terms which depend on the same power of ϵ, then, the E 0 and E 1 can be expressed as The keplerian elements can be transformed by the equinoctial elements, therefore, the variation of the keplerian elements due to the 122 slow-push action can be obtained Furthermore, the total variation in the mean anomaly between deflected and undeflected orbit is Therefore, the deflection vector δr f , t d , t f ; ae, ∆ae after the slow-push strategy can be calculated by Eq. (11). wherer,θ,ĥ are column vectors that define the radial, transverse, out-of-plane reference frame: ξ,η,ζ T are column vectors that define the B-plane reference frame: The deflection distance of the deflection action is measured by the projection of the deflection vector δr on the B-plane as Once the mean value and standard deviation of asteroid orbital elements and magnitude are given at the deflection time t d , the 132 mean value µ ξζ and covariance matrix µ ξζ ofb can be approximated by aforementioned UT technique as the following steps: Theb 133 should be reformed as a function T of orbital elementsae and the magnitudeH.
The mean value and covariance matrix ofX are Calculate the sigma points X and their weights , where j = 1, 2, ..., n; n ∈ N is the dimension of the state vector; β should be 2 for Gaussian distributions; and (. . . ) j means the j th column vector of the matrix. Obtain the set of the transformed sigma points Y (i) : 138 The mean value and the covariance are by using the weights and transformed sigma points as follows: Therefore, the Probability of Collision P ′ c after the deflection can be expressed as where P ′ c is computed by integrating the uncertainty ellipsoid, centered on the asteroid's mass point and projected on the B-plane, 142 over the closed region B ((0, 0) , R) defined by Earth's radius (assuming R E = 6378km in this paper). Patera's method Patera (2001) 143 is used for calculating P ′ c due to the fact that computational efficiency is particularly advantageous when large numbers of P ′ c 144 evaluations are performed.
with a bpa(e i j ) ∈ [0, 1] associated to each interval. Then the uncertainty set U is given by the Cartesian product U = I 1 × I 2 × ...I n 162 and we can define a focal element γ q = e 1J q (1) × e 2J q (2) × ...e iJ q (i) × ...e nJ q (n) with associated bpa(γ q ) = i bpa iJ q (i) where the vector J q has n components and contains a permutation of indexes j. We can now define the set A ν as: and the cumulative Belief and Plausibility associated to proposition in Eq.( 34): Once the Belief in a given value of ν is computed, the following multi-objective optimisation problem can be formulated in order to 167 maximise the Belief in the minimum achievable P ′ c : The optimal design vector and thresholds that yield a Bel = 1 for all possible u ∈ U can be computed by solving the following classic   the NEO, its mass and the associated uncertainties (see Figure 1). For each virtual impact scenario and deflection technology an  This section explains how we generated the data-set to train the machine learning model. The procedure involves three main 209 steps: the first step is to generate a variety of Virtual Impactors (VIs) with different close encounter geometries in order to cover as 210 many possible future scenarios as possible, the second step is to apply uncertainties to virtual impactors to form the virtual impact 211 scenarios, and the third step is to perform the robust optimisation for each virtual impact scenarios and label the resulting probability 212 of collision.

213
Step1 (generate virtual impactors): we assume the Earth orbit is circular, therefore, two necessary but not sufficient conditions on 214 the semi-major axis a, eccentricity e for virtual impactors are Fixing the semi-major axis a, eccentricity e and inclination i with Near Earth Object's actual value from the JPL Small-Body 216 Database, one independent element remaining to fix is the longitude of the ascending node node Ω of the asteroid's orbital plane   Table 2. For each VI, the warning time are randomly collected from [5, 10] years and uncertainty intervals e i j 227 are the random subsets collected from the uncertainty intervals Λ i j listed in Table 2, that is 228 {e i j |e i j ⊂ Λ i j , i = 1, 2, ..., 7; j = 1, 2} Finally, a total of 5,000 virtual impact scenarios are generated.

229
Step3 (perform robust optimisation and label the outcome): For each virtual impact scenario, we perform a robust optimisation, where d * is the robust solution of problem (40). Figure 5 shows the worst deflection effectiveness P ′ c(max) of five different deflection 233 strategies within 10 years warning time, where the x, y, z axis indicate the semi-major axis, eccentricity and diameter of the VIs 234 respectively. The P ′ c(max) will be used to determine whether the execution of deflection is successful.

235
In this paper, a 'successful deflection' is defined as a deflection such that the worst probability of collision P ′ c(max) is below 10 −2 , 236 that is: In other words, once the P ′ c(max) of a certain deflection strategy is less than 10 −2 , the deflection strategy will be labeled as 'successful', 238 otherwise it will be labeled as 'un-successful'. Figure 6 shows the percentage of the scenarios (among 5000 virtual impact scenarios)  ', 50, 100); the minimum number of samples required to be at a leaf node is (1, 10 −4 , 10 −7 ); the minimum 251 number of samples to split a node is (2, 20); the number of features to consider when looking for the best split is ('auto',0.5,'log2').

252
For each classifier, we use 80% of the samples for training and 20% for testing. The mean value of F1-score is employed to assess The points that tend to be blue in color indicate the lower P ′ c(max) , to some extent, symbolize the virtual impact scenarios that can be successfully deflected. Fig. 6: Percentage of the scenarios that can be deflected successfully. With the warning time less than 10 years, NED is the most powerful strategy that can defend over 80% asteroid impact scenarios, while the GT is the weakest strategy that can only defend 2.6% asteroid impact scenarios. Classifier, this is caused by the imbalance in the current data set. For example, the percentage of the scenarios labeled as 'successful' 263 for GT and IBS are less than 20%. In this case, the recall and precision scores of GT Classifier and IBS Classifier indicate that the 264 'successful' scenarios for GT and IBS tend to be mislabeled as 'un-successful'. Such imbalanced data set is intrinsic to the deflection 265 scenarios with a NEO diameter less than 1000 m and a warning time less than 10 years. An extended data set, which considers larger 266 NEOs and longer warning times, is needed to improve the performance of the IDSS.

267
The statistical analysis of the output of the current classifier also implies that the NED technology is reliably classified when 268 successful, however, since there is a small proportion of outputs in which the NED option is classified as unsuccessful even a small 269 percentage of false positives for the label "successful" drives the recall down to 44.19%. A symmetric consideration applies to the 270 GT and IBS because the percentage of times the answer is "successful" is small compared to the false positive labels "unsuccessful".

271
In both cases, however, the total number of correct answers is very high, thus, even with an unbalanced data-set, the IDSS returns the from Λ ·1 and Λ ·2 listed in Table 2, while two sources are equally reliable: The information and simulation results of the three cases are summarised in Table 4, where Case 1 indicates the Virtual Impactor 277 (VI) with the larger size and easier orbital accessibility, Case 2 indicates the VI with the smaller size and easier orbital accessibility IDSS is used to predict the suggested strategies, the total computational time for each case can be reduced to less than 2s, while the 282 prediction results are the same as the results from robust optimisation. Table 4 shows a reduction of nearly three orders of magnitude 283 in the computational time to assess the effectiveness of a group of deflection technologies. This reduction is proportional to the total 284 mission duration and is expected to further increase as more uncertainty and higher fidelity models are considered in the prediction.

285
For these three cases, the corresponding belief curves of five deflection strategies are shown in Figure 7a Table 2. Taking one focal element that includes the σ a collected from Source2 as an example, Figure 8 gives the 3-σ ellipses 289 (on the B-plane) corresponding to the best case P ′ c(min) and the worst case P ′ c(max) . The Case 1's VI has a diameter of 97m, which is  shows that NED and KI outperform slow-push strategies, at the same time, a slight advantage of LA is shown among the slow-push 294 strategies due to the fact that LA does not require to carry as much propellant on board, as IBS and GT, to deflect the NEO.

295
Note that a NED for small asteroid might not be a recommended solution for other reasons than the pure deflection efficiency.

296
In this paper we did not include additional constraints and considerations which would be required to select the nuclear option vs.   Thus the IDSS can be used to provide a rapid screening of all deflection options at an early stage of the decision making process 308 and multiple fast re-assessment of a given technology once more information is available. The IDSS can also be used to assess 309 the applicability of different technologies to a catalogue of NEOs and to quickly repeat the assessment when the impact risk in the 310 catalogue is updated.

311
The results in this paper are preliminary but demonstrate the potentiality of the combination of Dempster-Shafer theory of evidence 312 and machine learning to make informed decisions under mixed uncertainty on the most suitable deflection options in the case of 313 a threatening scenario. This work included only uncertainty on the mass of the asteroid and its ephemerides. A more complete 314 treatment would require the inclusion of uncertainty on the deflection action and on other physical properties of the asteroid. This is a critical aspect of the effectiveness and applicability of all deflection methods. The inclusion of uncertainty on the deflection action 316 and its sensitivity to the properties of the asteroid would change the classification in favour of more reliable methods.

317
Furthermore the only criterion used to classify is the probability of collision. This would need to be complemented by possible 318 operational or implementation constraints and more general considerations on technology readiness and viability. Finally the data-set 319 used to train the classifiers contains a lot more information on the mission profile, trajectory, transfer manoeuvres, launch date, 320 etc. In the current version of the IDSS this information is not returned but can be used to train a further model associated to the 321 classification of the technology so that an optimal mission is returned together with the associated technology. All these extensions 322 will be the subject of future work.