## 4 Numerical Simulations
### High-Earth Orbit Case  
First, consider a spacecraft orbiting on a high-Earth orbit (HEO). Assume that the motion of the spacecraft is confined to the equatorial plane, which allows the position \(r\) to be described by two scalar values \(r_{x}\) and \(r_{y}\), and the velocity \(v\) to be described by two scalar values \(v_{x}\) and \(v_{y}\). Thus, the state vector and equations of motion that describe the nonlinear dynamical system are written as follows [42]:  
\[ \begin{array}{l}\mathbf{x}=\left[\mathbf{r}^{T},\mathbf{v}^{T}\right]^{T}=\left[r_{x},r_ {y},v_{x},v_{y}\right]^{T}\\ \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})=\left[v_{x},v_{y},-\mu r_{x}r^{-3},-\mu r_{y}r^{-3} \right]^{T}\end{array}  \tag{27}\]  
where \(\mu=398,600\) is the gravitational constant of the central body (i.e., Earth in this paper). Take notice that although dynamics (27) are of the analytical form, it is still considered a black-box model by invoking the kinetic integral \(\mathbf{F}(t)=\int_{0}^{t}\mathbf{f}(\mathbf{x})d\mathbf{x}\). The elements of the nominal orbit are listed in Table 4, given by a semimajor axis of 35,000 km and an eccentricity of 0.2. Figure 2: Active sampling for one-dimensional numerical benchmark. \begin{table}
\begin{tabular}{c c c c c} \hline \hline _a/_km** & \(e\) & _i/_deg** & \(\mathbf{\Omega}\)/_deg_ & \(\mathbf{\omega}\)/_deg_ & _n/_deg_ \\ \hline
6596 & 0 & 0 & 0 & 0 & 0 \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Elements of the nominal orbit of LEO case. The initial state distribution is taken to be Gaussian, with the initial mean represented by the nominal orbit elements. Additionally, the initial covariance is set to be diagonal with standard deviations in the position elements of 1 km and the velocity elements of 1 m/s:  
\[\mathbf{P}_{0}=diag\Big{(}\Big{[}1\;\text{km}^{2},1\;\text{km}^{2},10^{-6}\;\text{ km}^{2}/\text{s}^{2},10^{-6}\;\text{km}^{2}/\text{s}^{2}\Big{]}\Big{)} \tag{28}\]  
To characterize the performance of the proposed AGMM-VSG method, both the AGMM-VSG and fifth-order CR methods are applied to the problem of predicting the evolution of uncertainty. The propagation time is one orbit period (around 65,164.833 s). The values of parameters to be used in AGMM-VSG are set as \(t_{f}=65\),164.833 s, \(dt=t_{f}/100\). The criterion for splitting the Gaussian components are \(v>10^{-4}\), where \(v\) is the nonlinearity index (Equation (15)) of the corresponding component. Parameters for active sampling are given as: \(n_{c}=200\), \(n_{r}=100\), \(n_{i}=50\), and \(n_{add}=10\). In addition, a Monte Carlo simulation using 1000 samples drawn from the initial distribution is performed. The position and velocity marginal PDF contours using AGMM-VSG and fifth-order CR are illustrated in Figures 3 and 4, respectively. The outputs of the Monte Carlo simulations are shown represented with a grey x. Depart from the scatter plot, the PDFs of each state variable are also illustrated. The PDFs of each state variable of AGMM-VSG and fifth-order CR are obtained by calculating the integral of the multivariate PDFs, while PDFs of each state variable of 1000 Monte Carlo samples are obtained using MATLAB function _ksdensity_. It can be seen from Figures 3 and 4 that, compared with fifth-order CR, the proposed method is better at representing the curvature exhibited by the Monte Carlo samples. For AGMM-VSG, only 65 Gaussian components and 100 real samples are requested. The computation cost of AGMM-VSG is a little higher than that of fifth-order CR, which uses a total of 33 evaluations (\(2n^{2}+1\) Cubature-points and \(n=4\)), while significantly outperforming fifth-order CR in terms of accuracy. Moreover, if use real samples for propagating, a total of 2145 real samples (65 components and 33 evaluations for each component) are needed, which indicates the signification of the utilization of virtual samples for reducing computational expense. agreement between two distributions.