A method for direct in-space thrust estimation from low-acceleration orbital maneuvers

This paper presents a method for performing direct in-space thrust estimation for low-acceleration propulsion systems using measurements of the spacecraft’s position taken during an orbital maneuver. The method is based on the ensemble Kalman update which does not require linearization of the spacecraft dynamics nor does it require Gaussian distributions for parameter uncertainties and measurement noise, allowing for a more general approach to thrust estimation. In addition, modeling error, such as that caused by the truncation of a spherical-harmonics representation of the Earth’s gravitational field, can be explicitly accounted for by representing the error with Gaussian processes. Simulated experiments show that uncertainty in the propulsive acceleration magnitude on the order of 0.1 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}m/s2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} (1σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document}) can be achieved at an orbit altitude of approximately 410 km with temporally-sparse measurements even in the presence of uncertain atmospheric drag, and Monte Carlo analysis demonstrates the consistency of the inference results. Trends in the estimate of the propulsive acceleration with the true acceleration value are explored empirically and theoretically in order to allow for generalization of the results. The outcome of this work is a systematic approach to direct in-space thrust estimation that can support the final steps of development for future in-space electric propulsion systems or the calibration of a thruster during a mission.


Introduction
During the in-space demonstration of a new propulsion system or calibration of the propulsion system performance during a mission, a direct measurement of the thrust produced by the propulsion system is usually desired such that the performance in space can be compared to the performance measured in a laboratory environment or predicted by numerical models.Direct in-space thrust measurements typically rely on using the propulsion system to induce an angular or linear acceleration on the parent spacecraft and measuring the resulting dynamic response.Both of these approaches were demonstrated in the Space Electric Rocket Test (SERT) missions: SERT-1 configured the propulsion system to produce torques on the spacecraft and measured the change in the spacecraft's spin rate [1] while SERT-2 configured the propulsion system to produce linear accelerations on the spacecraft and measured the change in the spacecraft's orbital radius [2].Both spacecraft also carried on-board accelerometers in order to take instantaneous measurements of the acceleration produced by their respective propulsion systems.
In many cases, measurement of the thrust through linear accelerations and orbital maneuvers may be the desired approach as it imposes the least amount of requirements on the propulsion system architecture.For missions that are not dedicated technology demonstration missions, measuring thrust through angular accelerations by offsetting the thrust vector from the center of mass of the spacecraft requires gimbaling mechanisms or that multiple thrusters be used throughout the primary mission.For some spacecraft this is possible with the desired propulsion system architecture, as was the case with the Mars Cube One spacecraft [3], but is not generally true.Direct measurement of the propulsive acceleration through an on-board accelerometer that can provide accurate and precise measurements of the propulsive acceleration takes away from valuable payload mass and volume, particular for the low accelerations (order 10 µm/s 2 ) expected when using an electric propulsion system.
These problems are exacerbated for electric micropropulsion systems used onboard small spacecraft which feature reduced budgets and operational support throughout the mission.Development of electric micropropulsion systems is also an extremely active research area [4]; many low-thrust propulsion systems for small spacecraft such as ionicliquid electrospray thrusters [5,6], film-evaporating water microcapillaries [7], miniaturized ion thrusters [8], and plasma thrusters [9,10] are actively under development and will likely have in-space demonstrations in the near future.Since these propulsion systems are designed for small spacecraft, it is likely that the in-space demonstration will be conducted on a dedicated technology demonstration mission with a CubeSat.
The combination of low-thrust propulsion with small spacecraft such as CubeSats creates a unique challenge for in-space thrust measurements.The low cost and size of the spacecraft likely prohibits the use of accelerometers with precision high enough to directly measure the thrust.Even with a dedicated technology demonstration mission where the propulsion architecture is less constrained, the small size of the spacecraft can also present issues for measurements relying on angular accelerations due to uncertainty in the lever arm of the propulsion system.For a CubeSat, unless the propulsion system is mounted on an extendable boom, the maximum lever arm will be on the order of a few centimeters.Therefore, uncertainty in the lever arm of the propulsion system of only a few millimeters will lead to errors in the thrust measurement on the order of 10%.Uncertainty in the lever arm can come from many sources including the thruster itself.For example, spatial variation of emission across an ionic-liquid electrospray thruster [11][12][13] could lead to variations of the center of thrust of the thruster on the order of a few millimeters, placing an inherent limit on the achievable precision of the thrust estimate.
Due to these limitations, while methods for estimating propulsion system thrust from orbital maneuvers may be desirable in general, these methods may be required for direct in-space thrust estimation of electric micropropulsion systems.One of the simplest approaches for estimating thrust from an orbital maneuver is to perform a circular orbit raise.An analytical approximation for the radial position of the spacecraft can be derived where r is the radial position of the spacecraft, v is the spacecraft's speed, a p is the pro- pulsive acceleration, and t is time [14].By measuring the spacecraft's radial position over time and inverting Eq. 1, an estimate of the propulsive acceleration can be determined, from which an estimate of the thrust can be obtained.This approach was used during the SERT-2 mission [15].However, there are many limitations for measuring thrust in this manner.First, it requires that the propulsion system be active for the entire duration of the maneuver which may not be possible, especially for power-limited spacecraft such as CubeSats.Second, only a point-estimate of the thrust is obtained and the uncertainty in the estimate can be difficult to quantify [16], especially in low-Earth orbit where the analytical approximation is a poor representation of the spacecraft radius due to oscillations primarily from Earth's oblateness.Lastly, it is restricted to maneuvers where the spacecraft starts in a near-circular orbit and cannot be generalized to other maneuvers.These limitations mean that changes to the spacecraft's radial position can be used to check that the thrust is close to what is expected [17,18], but an alternative method for estimating thrust from orbital maneuvers is desired.
Augmented Kalman filters have been considered for maneuver detection of uncooperative satellites [19,20] and a linear batch filter has been evaluated for estimating thrust from a low-thrust maneuver [21,22].However, both of these approaches require linearizing the spacecraft dynamics in order to propagate the spacecraft's state over time.Due to the significant nonlinearity of orbital dynamics, a linearization of the spacecraft's dynamics requires that measurements be taken frequently in order to produce reliable results.Measurement intervals of 1 s [19] and 30 s [22] are used, which may be difficult to achieve on resource-constrained spacecraft such as CubeSats while also powering an electric propulsion system.In addition, the impact of modeling error during the inference step is rarely accounted for.While there have been some attempts to account for the modeling error of Earth's gravitational field, inadequate accounting for this modeling error can introduce biases in the estimate that exceed the reported uncertainty bounds of the inference method [22].In this case, the reported uncertainty of the inference method is unreliable, and does not capture the true uncertainty in the thrust estimate.
In this work, an inference approach is considered where the thrust is inferred based on measurements of the spacecraft's position taken over an entire maneuver.Specifically, for a system that can be parameterized by a parameter vector ψ and produces a meas- urement vector , Bayes' rule gives that the posterior probability density of the parameter vector given the measurement vector is proportional to the product of the likelihood of observing the measurement vector given the parameter vector with the prior knowledge of the parameter vector The goal of this inference method is to generate samples from the posterior distribution, P(ψ| ) .By using a sampling-based approach the likelihood probability, P( |ψ) , and prior probability, P(ψ) , can be any distribution rather than being restricted to Gauss- ian distributions.Additionally, by using an inference approach, measurements over the (1) (2) P(ψ| ) ∝ P( |ψ)P(ψ) entire trajectory are considered simultaneously, allowing the nonlinearities of the spacecraft's dynamics to be explicitly accounted for.Although this method is applicable to any level and duration of thrust, the focus of this work will be on cases with low acceleration as they represent the more-difficult scenario.
This paper develops a method for solving the above inference problem, more explicitly defined in Section "Problem statement".The method is based on the ensemble Kalman update with modeling error incorporated through process noise represented with Gaussian processes.The choice of algorithm is intentional in its inherent simplicity as well as flexibility to different measurement times, maneuvers, and other parameters of the inference problem.Sections "Problem statement"-"Modeling error" develop different aspects of the method while Sections "Implementation"-"Example" apply the method to an example thrust inference problem.The intention of this paper is to outline a systematic approach to direct in-space thrust estimation that can be used for future inspace performance characterization of electric propulsion systems, particularly electric micropropulsion systems.

Problem statement
Given m measurements of the spacecraft's position at times t ∈ {t 1 , ..., t m } , the measure- ment vector, , is a 3m × 1 vector containing the Cartesian x, y, and z, positions of the spacecraft at each measurement location The problem addressed in this work, is given an observed measurement vector, * , infer the initial Cartesian position and velocity of the spacecraft along with the average atmospheric drag acceleration magnitude, a d , and the average propulsion system acceleration magnitude, a p , averaged over the entire maneuver period.Only the magnitudes for the atmospheric drag and propulsive accelerations are considered here, as their directions are assumed to be known.These parameters form an 8 × 1 parameter vector The parameter and measurement vectors are related through a forward model where f (ψ) is a function that takes the parameter vector, propagates the spacecraft's position and velocity forward in time, and returns the spacecraft's position at the desired measurement times.In this work, both the measurement and parameter vector are assumed to be expressed in an Earth-centered, Earth-fixed reference frame, but the same method can be applied to alternative reference frames.ν is a 3m × 1 vector representing the measurement noise where each of the ξ i are independent, normally-distributed random variables with mean of zero and standard deviation of σ : ξ i ∼ N (0, σ 2 ).
(3) = x(t 1 ) ... x(t m ) y(t 1 ) ... y(t m ) z(t 1 ) ... z(t m ) T By performing the inference, the average propulsion system acceleration in the presence of uncertain atmospheric drag can be determined which allows for a characterization of the propulsion system performance.Note that the inferred propulsive and atmospheric drag accelerations are assumed to be spatially and temporally constantthis restriction will be further discussed in Section "Dynamics models".Alternative parameter vectors that allow for spatial and temporal variation in the propulsive and atmospheric drag accelerations can be used (e.g.inferring the coefficients of a polynomial function that describes the temporal change in the propulsive acceleration).However, in practice it has been found that these variations are difficult to discern relative to inferring a single average value due to the small magnitude, and therefore small dynamic influence, of the propulsive and atmospheric drag accelerations.

Ensemble Kalman update
The ensemble Kalman Update [23] is a form of approximate Bayesian inference meaning that the distribution of samples it generates approximates the true posterior distribution.The idea behind the ensemble Kalman update is the same as the update step of a standard Kalman filter where parameters of inference are updated based on a linear update rule where K is the optimal Kalman gain, ψ + is the posterior estimate for the parameter vector, (ψ − ) is the simulated measurement vector based on the prior estimate for the parameter vector ( ψ − ), and * is the observed measurement vector.For linear-Gauss- ian models, the Kalman gain can be calculated analytically based on the measurement model and prior covariance.However, for models that contain nonlinearities or non-Gaussian distributions, an analytical form of the optimal Kalman gain is generally not possible to obtain.
Instead, the ensemble Kalman update attempts to approximate the optimal Kalman gain through Monte Carlo methods.The key idea is that in the linear-Gaussian case, the optimal Kalman gain is calculated from where � ψ, is the covariance between the parameter and the measurement vectors, and is the variance of the measurement vector.By generating an ensemble of initial parameter vectors, based on the assumed prior distribution of the inference parameters, the measurement vector for each ensemble member can be simulated from the forward model of the system, Eq. 5.With a collection of parameter vectors and corresponding measurement vectors, � ψ, and can be numerically estimated which allows for an estimate of the optimal Kalman gain to be obtained.
Algorithm 1 shows the ensemble Kalman update routine.The ensemble, , consists of a collection of parameter vectors.For a given number of ensemble members, the parameter vector for each member of the ensemble is generated by randomly sampling a parameter vector from the prior distribution, P(ψ) .Then, the forward model of the system, f (ψ) , is used to generate a simulated measurement vector given the parameter (7) vector for each ensemble member along with additive measurement noise, ν .Given the collection of parameter vectors and simulated measurement vectors, the covariances � ψ, and can be numerically estimated which allows for a numerical estimate of the optimal Kalman gain to be obtained.With the estimate of the optimal Kalman gain, the parameter vector of each ensemble member can be updated using the linear update rule in Eq. 7. The posterior distribution of the parameter vectors in the ensemble is then an approximation of the true posterior distribution for the inference parameters.

Dynamics models
To actually use the ensemble Kalman update, the forward model needs to be defined.In this work, the forward model propagates the spacecraft's position and velocity forward in time based on the initial spacecraft state and average propulsive and atmospheric drag accelerations, and returns the spacecraft's position at the desired measurement times.The key aspect of the forward model is therefore the underlying dynamics model.Two models of the dynamics in Earth orbit are used in this work: a high-fidelity model for simulating the observed data and a low-fidelity model used in the forward model for the ensemble Kalman update.The high-fidelity model is the General Mission Analysis Tool [24] and uses a 70th-order spherical-harmonics representation of the Earth's gravitational field with Joint Gravity Model 3 coefficients [25], the MSISE90 atmospheric model [26], solar radiation pressure, and point-mass representations of the Sun, Moon, and other planets.This model is only used to simulate the data that would be collected on-orbit and is not used in the ensemble Kalman update during the inference step.
The low-fidelity model only accounts for the gravitational field of the Earth and atmospheric drag, and does not consider solar radiation pressure or the influence of other planetary bodies as they are negligible relative to the gravitational and drag accelerations.In addition, the spherical-harmonics representation of the Earth's gravitational field is truncated to 30th order and the atmospheric drag acceleration is assumed to be spatially and temporally constant.This low-fidelity model is used in the forward model for the ensemble Kalman update, f (ψ) .In principle, the high-fidelity model can be used in the inference step.However, the use of a low-fidelity model offers two advantages: it reduces the computational complexity of the model and reduces the number of parameters that need to be inferred.In addition, the low-fidelity model alleviates potential issues with overfitting a specific model to the measurements that could bias the estimated propulsive acceleration, particularly for modeling the atmospheric drag.
With the low-fidelity model, the acceleration vector of the spacecraft in an Earth-centered, Earth-fixed frame becomes where r and v are the spacecraft's position and velocity vectors respectively, û is the unit vector of the propulsive acceleration vector, v is the unit vector of the spacecraft's veloc- ity vector, g is the local gravitational acceleration, and is the rotational frequency vec- tor of the Earth.It is assumed that the atmosphere rotates with the Earth, therefore the appropriate velocity direction for determining the atmospheric drag acceleration vector is the spacecraft's velocity in the Earth-centered, Earth-fixed frame.The final term of the acceleration, n(t) , is process noise added to the simulation in order to represent mod- eling error, and will be further discussed in Section "Modeling error".
Only eight parameters are required in order to simulate a spacecraft trajectory in an Earth-fixed, Earth-centered frame using the acceleration vector in Eq. 9: the initial position and velocity of the spacecraft, the average atmospheric drag acceleration, and the average propulsive acceleration.These parameters are therefore the parameters of inference and form the parameter vector in Eq. 4. The forward model for the ensemble Kalman update is then just a propagation of the spacecraft's position and velocity over time for a particular parameter vector given the acceleration vector of Eq. 9. Measurements are obtained by taking the spacecraft's position at the desired measurement times and adding additive measurement noise, per Eq. 5, in order to form the measurement vector in Eq. 3. Noise in the measurements is assumed to be added to all three axes in Cartesian coordinates with the noise along each axis independent of the other axes.This noise model was chosen for simplicity, but other noise models, including those that contain nonlinearities or correlation between noise in the different axes, can be used.

Modeling error
One of the key aspects in inferring the probability distribution for the propulsive acceleration is the handling of modeling error.In this problem, modeling error will be caused by errors in the acceleration applied to the spacecraft from different sources, and will determine the process noise, n(t) , in Eq. 9.For the choice of low-fidelity model used in this work, the dominant source of modeling error will be caused by truncation of the spherical-harmonics representation of the Earth's gravitational field.Making the Earth's gravitational field the dominant source of modeling error is an intentional choice due to the availability of accurate high-order models for comparison.Since the Earth's gravitational field is the dominant source of modeling error, the process noise can be easily determined by comparing the low-fidelity model to a known higher-fidelity representation of the Earth's gravitational field.The resulting process noise will be conservative, but more accurately represent the error between ( 9) the low-fidelity and high-fidelity models relative to using a different source of modeling error as the dominant source.In particular, even high-fidelity models of atmospheric density, such as the MSISE90 atmospheric model [26], have significant (10%) and unknown error which creates a challenge in quantifying the corresponding process noise.Table 1 shows the approximate order of magnitude for different sources of modeling error in an Earth-fixed, Earth-centered reference frame for a near-circular orbit with altitude of 410 km.When calculating the acceleration errors due to atmospheric drag and solar radiation pressure, the spacecraft is assumed to have a drag reference area of 0.1 m 2 , coefficients of drag and reflectivity of 1.0, and mass of 4 kg, representative of a 3U CubeSat.
While the truncation of the spherical-harmonics representation of the Earth's gravitational field is intentionally chosen to be the dominant source of modeling error, other sources of error can still have an effect.In particular, for a near-circular orbit the atmospheric drag acceleration will predominantly act in the along-track direction.If the propulsive acceleration vector is also fixed to the along-track direction, as it would be for an orbit raising or lowering maneuver, errors in estimating the atmospheric drag acceleration can show up as a bias in the estimate of the propulsive acceleration.Since the atmospheric drag acceleration is assumed to be spatially and temporally constant in this work, slight errors in the atmospheric drag acceleration between the low-fidelity and high-fidelity dynamics models at any instant in time are guaranteed to exist.This bias is expected to be small as uncertainties in the along-track acceleration due to truncation of the spherical-harmonics representation of the Earth's gravitational field are designed to be the dominant error source, and Monte Carlo analyses in Section "Monte Carlo analysis" demonstrate that any potential biases from modeling of the atmospheric drag acceleration are small relative to the reported uncertainty from the ensemble Kalman update.
Figure 1 shows the acceleration error in the along-track direction caused by truncating the spherical-harmonics representation of the Earth's gravitational field from a 70th-order to a 30th-order spherical-harmonics representation for a 410 km altitude, near-circular orbit.The magnitude of the acceleration error has some cyclical behavior with period close to the 90-minute orbital period.However, the error qualitatively appears to be random noise.Furthermore, Fig. 2 shows a histogram of the acceleration error which shows that the magnitude of the acceleration error is near-Gaussian.These two aspects combined suggest that generating the process noise, n(t) , for Eq. 9 from a Gaussian process may be a reasonable representation of the modeling error due to truncating the spherical-harmonics representation of the Earth's gravitational field.
In order to generate realizations of the process noise from a Gaussian process, the covariance between samples in the process noise needs to be defined.Figure 3 shows the sample auto-correlation between samples of the along-track acceleration error along with an exponentiated-quadratic curve fit for the kernel of the Gaussian process where τ is the scale of the Gaussian process, t is the time between samples or lag, and L is the length of the Gaussian process.The exponentiated quadratic is able to represent the sample auto-correlation for a lag of approximately one minute and lower, but cannot (10) Fig. 1 Along-track gravitational acceleration error for a 410 km, near-circular orbit Fig. 2 Probability density of along-track gravitational acceleration error for a 410 km altitude, near-circular orbit.Dash-dot line represents a Gaussian fit to the histogram data represent the auto-correlation for longer lags.Instead, samples that are separated by one minute or greater are assumed to be uncorrelated.Future work will explore the use of different kernels for the Gaussian process in order to better represent the true sample auto-correlation.
With the scale of the kernel determined from the variance in samples from Fig. 2 and the length of the kernel determined from Fig. 3, realizations of the Gaussian process can be generated that provide a reasonable approximation of the expected modeling error.Note that this approach differs from prior approaches, such as that in Ref. [22], in that the correlation in the process noise is accounted for rather than assuming that the process noise is uncorrelated white noise.Since it can be seen from Figs. 1 and 3 that the modeling error is not uncorrelated white noise, generating the process noise from a Gaussian process will provide a more-accurate representation of the expected modeling error.Figure 4 shows one such realization for the alongtrack gravitational acceleration error.Qualitatively, it looks similar to the true error in Fig. 1, but with a greater content of high-frequency noise.The increase in highfrequency noise is caused by the poor representation of the sample auto-correlation for lags greater than one minute.
Although only the modeling error in the along-track direction is shown here, the same process can be applied to the model error of the gravitational acceleration in the radial and cross-track directions.A total of three Gaussian processes are defined, one for each axis, and serve to represent the process noise in Eq. 9. Since the Gaussian processes are defined in a local-vertical, local-horizontal frame but the acceleration of Eq. 9 is defined in a Cartesian frame, the process noise at any given point in time is generated in the local-vertical, local-horizontal frame and then rotated to the Cartesian frame based on the spacecraft's current position and velocity.The generated noise in the Cartesian frame is then combined to form the process noise vector, n(t) , in Eq. 9.

Implementation
For this work, a simple maneuver is considered where the spacecraft maintains its thrust axis aligned with the along-track direction where ĥ is the unit vector of the orbit's angular momentum, calculated in an Earth-cen- tered inertial frame, and r is the unit vector of the spacecraft's position.There are no restrictions in the method on the direction of the thrust vector.However, a thrust vector aligned with either the angular or velocity directions will produce a more-significant change to the spacecraft's position, due to the increase or decrease in orbit semi-major axis, relative to one aligned with the radial or cross-track directions.Larger changes in the spacecraft's position improve the observability of the propulsive acceleration as it will be easier to differentiate the changes due to propulsive acceleration from those due to measurement noise.The maneuver is also broken into two distinct phases of equal length: a powered phase where the propulsion system maintains a constant thrust and a decay phase where the propulsion system is off and the spacecraft's orbit is allowed to decay due to atmospheric drag.This maneuver is used to aid in differentiating the effects of atmospheric drag and propulsive thrust, and allow for inference of both the average atmospheric drag and propulsive accelerations.Measurements of the spacecraft's position are collected at a fixed interval of 10 min over the entire trajectory and used to infer to the initial state of the spacecraft, the average atmospheric drag acceleration, and the average propulsive acceleration.There is no requirement that the measurements have to be taken at a fixed time interval, and the effects of uneven sampling will be explored in future work.Gaussian measurement noise with 3σ uncertainty of 10 m is added to the measurements along each axis in order to represent measurements taken by a Global Positioning System antenna.
For all results, the true spacecraft trajectory is simulated from the initial orbital parameters given in Table 2, representative of deployment from the International Space Station.The initial spacecraft mass was set to 4 kg, the coefficients of drag and reflectivity (11) û = ĥ × r Fig. 4 Realization of a Gaussian process for representing the along-track gravitational acceleration error for a 410 km altitude, near-circular orbit were both set to 1.0, and the drag and solar radiation pressure reference areas were both set to 0.1 m 2 .For the MSISE90 atmospheric model, the F 10.7 solar flux and average solar flux were both set to 150 solar flux units while the geomagnetic index was set to 3. The propulsion system was assumed to provide a constant thrust with specific impulse of 1000 s.While propellant mass flow was simulated, due to the relatively high specific impulse and relatively short maneuver lengths, the impact of propellant mass flow on the trajectory is negligible.It is finally assumed that the thrust vector is perfectly aligned with the desired direction for the entire maneuver.The impact of uncertain thrust vector direction will be explored in future work, but is expected to have a negligible effect due to the relatively accurate pointing (<5 degrees) of small spacecraft attitude control systems-a constant 5 degree error in the thrust vector will only lead to a 0.4% decrease in the along-track thrust.
In order to determine the appropriate parameters for the Gaussian process used to generate realizations of the process noise, an orbit close to the expected orbit of the spacecraft is simulated.This orbit does not need to exactly match the actual orbit of the spacecraft, but should be qualitatively similar such that reasonable estimates of the modeling error can be determined.Here, the nearby orbit is determined based on the initial measurement of the spacecraft's position.Based on the simulation of the nearby orbit, the acceleration error in the radial, along-track, and cross-track directions due to truncating the spherical-harmonics representation of the Earth's gravitational field can be determined.The process outlined in Section "Modeling error" was used in order to determine the scale and length of the respective kernels.Table 3 shows the scale and length of the kernels in all three axis based on a simulation of a 410 km altitude, nearcircular orbit.This process only requires a single simulation with a high-fidelity model of the Earth's gravity field as opposed to the potentially thousands of simulations used in the ensemble Kalman update, thereby dramatically reducing the computational cost of performing the inference step.For the prior distributions of the initial spacecraft position, Gaussian distributions with mean given by the initial measurement of the spacecraft's position and 3σ uncer- tainty of 10 m are used.Since it is assumed here that no direct measurement of the spacecraft's velocity is available, the distributions for the initial spacecraft velocity are assumed to be uniform with mean determined from a best fit between the position of the spacecraft given by the first two measurement locations and under the assumption of no propulsive or atmospheric drag acceleration.The range of the uniform distributions was assumed to be 2 m/s, which empirically was found to best represent the uncertainty in the prior velocity estimate from the best-fit method.The prior distributions for both the atmospheric drag and propulsive accelerations are both assumed to be uniform and range from zero to maximum values of 7 µm/s 2 and 40 µm/s 2 respectively in all cases considered in this work.

Example
Figure 5 shows the simulated radial position of the spacecraft during a 16-hour maneuver using the high-fidelity model (GMAT) where the propulsion system is on for the first eight hours with a thrust of 100 µ N (acceleration of 25 µm/s 2 ) and off for the last eight hours.Measurements of the spacecraft's full three-axis position are taken once every 10 minutes for a total of 97 measurement times throughout the entire maneuver.Figure 6 shows the resulting marginal posterior distribution of the propulsive acceleration using the ensemble Kalman update approach with an ensemble size of 2,500.The posterior distribution is near-Gaussian with mean of 24.809 µm/s 2 and standard deviation of 0.277 µm/s 2 giving a 0.8% error in the mean value and 3.3% 3σ uncertainty.
In addition to an estimate of the propulsive acceleration, the ensemble Kalman update provides an estimate of the average atmospheric drag acceleration throughout the trajectory.Figure 7 shows the marginal posterior distribution of the atmospheric drag acceleration.The true average atmospheric drag acceleration throughout the trajectory was approximately 3.520 µm/s 2 .The mean of the inferred distribution from the ensemble Kalman update is 3.348 µm/s 2 and the standard distribution is 0.196 µm/s 2 giving a 4.9% error in the mean value and 16.7% 3 σ uncertainty.Figure 8 shows a scatter plot of the propulsive and atmospheric drag accelerations for the posterior ensemble members.There is a clear positive correlation between the two accelerations, which is to be expected; trajectories with a larger atmospheric drag acceleration require a larger propulsive acceleration in order to explain the observed spacecraft position.

Monte Carlo analysis
A Monte Carlo analysis was conducted in order to check the consistency of the ensemble Kalman update approach.In the Monte Carlo analysis, simulated data was generated using the high-fidelity model with the maneuver time randomized between 2-16 hours (thruster firing times of 1-8 hours) and propulsive acceleration randomized between 0-40 µm/s 2 for each sample.Figure 9 shows the distribution of the error in the mean value of the inferred propulsive acceleration normalized by the standard deviation from the ensemble Kalman update from 500 Monte Carlo samples.For all samples, the spacecraft was initialized with the initial orbital parameters in Table 2.If the ensemble Kalman update were exact, the distribution should be a standard normal distribution (mean of zero, standard deviation of one).The actual distribution of the normalized mean propulsive acceleration errors has a mean of 0.177 and standard deviation of 0.630.The estimate of the propulsive acceleration from the ensemble Kalman update therefore has a slight bias, but the reported uncertainty is conservative.
The discrepancy in the mean and distribution of the normalized propulsive acceleration error in the Monte Carlo analysis is hypothesized to be due to the fact that the actual model error in the Earth gravity field for a fixed initial spacecraft orbital parameters is not true process noise as it is considered in Section "Modeling error".To test this hypothesis a second 500-sample Monte Carlo analysis was performed with the same Fig. 8 Scatter plot of propulsive versus atmospheric drag acceleration for the posterior ensemble Fig. 9 Distribution of the error in the mean propulsive acceleration estimate normalized by the inferred standard deviation from the ensemble Kalman update.In this case, the initial orbital parameters of the spacecraft were held constant.Dash-dot line represents a standard normal distribution parameters as before, but with the right-ascension of the ascending node additionally randomized from 0-360 degrees.Figure 10 shows the distribution of the error in the mean value of the inferred propulsive acceleration normalized by the standard deviation from the ensemble Kalman update with the randomized right-ascension of the ascending node.The mean offset is -0.115 and the standard deviation is 0.722 demonstrating slightly better accuracy of the reported mean and standard deviation for the estimated propulsive acceleration from the ensemble Kalman update.The remaining error may be caused by the discrepancy in representing the sample autocorrelation in Fig. 3 and will be explored in future work.However, this Monte Carlo analysis demonstrates the ability to generalize the calibrated Gaussian process parameters to other orbits, in this case those with different right-ascensions of the ascending node.Future work will explore the limits for generalizing the Gaussian process parameters calibrated for particular orbits.
It is worth noting that this Monte Carlo analysis demonstrates that the inferred distribution for the propulsive acceleration from the ensemble Kalman update is slightly conservative in that the reported standard deviation is larger than the actual standard deviation of the error between the mean and true propulsive accelerations.This is in contrast to prior approaches to accounting for modeling error in methods for direct in-space thrust estimation that lead to overly-optimistic predictions where the true propulsive acceleration does not fall within the predicted posterior distribution [22].With the ensemble Kalman update approach and accounting for modeling error with a Gaussian process, this Monte Carlo analysis shows that the propulsive acceleration is safely within three standard deviations of the reported mean value and that there is only a small ( ∼0.1σ ) bias in the estimator.
In addition to the Monte Carlo analysis demonstrating the extensibility of the Gaussian-process approach in approximating the modeling error to different orbits, trends in the posterior uncertainty in the propulsive acceleration can be explored versus maneuver parameters in order to generalize inference results to other maneuvers.Figure 11 shows the standard deviation of the marginal posterior distribution of the propulsive Fig. 10 Distribution of the error in the mean propulsive acceleration estimate normalized by the inferred standard deviation from the ensemble Kalman update.In this case, the initial right-ascension of the ascending node was randomized.Dash-dot line represents a standard normal distribution acceleration versus different values of the true propulsive acceleration.In all cases, the maneuver was fixed to eight hours where the propulsion system is active followed by eight hours of decay due to atmospheric drag with measurements taken every 10 minutes.Within numerical noise, the standard deviation is constant with propulsive acceleration.
This trend can be explained based on recent developments that analytically approximate the sensitivity of the spacecraft's position with respect to the applied propulsive acceleration, and show that the sensitivity is independent of the propulsive acceleration [27].Therefore, the observed variance in spacecraft position does not contain any information about the actual acceleration.Since the variance of the posterior distribution for the propulsive acceleration is determined by the measurement variance, it is unsurprising that the posterior distribution of the propulsive acceleration is independent of the true acceleration value since the measurement variance is also independent of the true acceleration value.Additional trends of the posterior uncertainty in the propulsive acceleration versus other maneuver parameters, such as maneuver duration and measurement interval, can be explored, but are heavily dependent on the process noise used in the inference process.As such, results regarding these trends are not presented here, and will need to be evaluated on a case-by-case basis.

Conclusion
This work shows that inference of the propulsive acceleration for low-thrust propulsion systems based on an ensemble Kalman update and measurements of the spacecraft's position over time can provide an accurate and precise estimate of the true propulsive acceleration in the presences of uncertain drag forces and modeling error.For a maneuver composed of eight hours where the propulsion system is active followed by eight hours of decay due to atmospheric drag at an altitude of approximately 410 km, an estimate of the propulsive acceleration with uncertainty of 0.831 µm/s 2 (3σ ) can be obtained when measurements are taken very 10 minutes.For a 4 kg spacecraft, such as a 3U CubeSat, this corresponds to 3.32 µ N of uncertainty in the thrust while for a 20 kg spacecraft this corresponds to 16.6 µ N of uncertainty.In addition, both empirical and theoretical exploration of the trend in the uncertainty in the propulsive acceleration versus true propulsive acceleration found that the absolute level of posterior uncertainty in propulsive acceleration is independent of the true propulsive acceleration.The method developed in this work has recently been applied to in-space performance characterization of an iodine electric micropropulsion system [28].
The maneuver used for the example in this work considers a single continuous thruster firing followed by an equal amount of time of decay due to atmospheric drag.However, the use of an ensemble Kalman update for thrust inference is not limited to this maneuver.Maneuvers where the spacecraft has multiple, separated segments where the propulsion system is active can be used.The only requirements are that the timing of when the propulsion system was turned on and off needs to be known such that the trajectory can be properly simulated during inference, and that there be an approach to differentiating the propulsive and atmospheric drag accelerations.Differentiating the propulsive and atmospheric drag accelerations can be accomplished with distinct phases of the maneuver where the orbit is allowed to decay due to atmospheric drag, as is done in here, or by aligning the thrust vector of the propulsion system with the radial or cross-track directions as is considered in [22].Considerations for the optimal thrust vector direction and maneuver design will be explored in future work.
The primary source of uncertainty in the propulsive acceleration estimate is due to modeling error.The impact of modeling error as it is incorporated in this work depends on the fidelity of the model used during inference as well as a subset of the orbital parameters.In this work, the fidelity of the model was intentionally selected such that the error in the gravitational acceleration was the dominant source of error, and can be accounted for through representing the acceleration error as a Gaussian process.It is of vital importance that the modeling error be considered and accounted for properly.Without considering model error and accounting for it during inference, it is likely that the resulting propulsive acceleration estimate will be overly optimistic where the true propulsive acceleration will not fall within the predicted posterior distribution.The Monte Carlo analysis of this work demonstrates that despite the relative simplicity of the dynamics model used during inference, the resulting posterior distributions for the propulsive acceleration are slightly conservative, and the true propulsive acceleration is safely within three standard deviations of the reported mean value.

Fig. 3
Fig. 3 Sample auto-correlation of along-track gravitational acceleration error for a 410 altitude, near-circular orbit.Solid line represents an exponentiated-quadratic fit to the data

Fig. 5
Fig. 5 Radial position of the spacecraft throughout the example maneuver.Markers indicate measurement locations

Fig. 6 Fig. 7
Fig. 6 Marginal posterior distribution of the propulsive acceleration obtained using the ensemble Kalman update.Dash-dot line represents a Gaussian fit to the histogram data

Fig. 11
Fig. 11 Standard deviation of the posterior marginal distribution of the propulsive acceleration versus true acceleration value

Table 1
Approximate order of magnitude for different sources of modeling error in an Earth-fixed, Earth-centered reference frame

Table 2
Initial spacecraft orbital parameters

Table 3
Gaussian process parameters for representing acceleration modeling error along each axis