A differential algebra technique to handle observation errors

The orbit determination is based on the observations of celestial bodies, which are affected by uncertainties due to instruments precision. Differential algebra can be used as a non-expensive computational tool to handle the observation noise. In the following, differential algebra is introduced. Then, its application in the orbit determination field is briefly explained.

Differential Algebra

Let us consider set of ordered pairs $(q_0,q_1)$, such that $q_0,q_1$ and $t$ are real numbers. Define addition, scalar multiplication and vector multiplication as follows:

\[\begin{gather} (q_0,q_1)+(r_0,r_1)=(q_0+r_0,q_1+r_1)\\ t\cdot(q_0,q_1) = (tq_0,tq_1)\\ (q_0,q_1)\cdot(r_0,r_1)=(q_0r_0,q_0r_1+q_1r_0) \end{gather}\]

being $t$ a real scalar. The ordered pairs with the introduced arithmetic are called $_1D_1$ and the three operations form an algebra. Now, let introduce a map $ \partial$ from $_1D_1$ into itself, defined by:

\[\partial(q_0,q_1) = (0,q_1)\]

Note that:

\[\begin{gather} \partial((q_0,q_1)+(r_0,r_1))=(0,q_1+r_1)=(0,q_1)+(0,r_1)= \partial(q_0,q_1)+ \partial(r_0,r_1) \\ \partial((q_0,q_1)\cdot(r_0,r_1)) =(0,q_0\cdot r_1+q_1\cdot r_0)=(0,q_1)\cdot(r_0,r_1)+(q_0,q_1)\cdot(0,r_1)=\\ = \partial(q_0,q_1)\cdot(r_0,r_1)+(q_0,q_1)\cdot \partial(r_0,r_1)\end{gather}\]

This holds for all $(q_0,q_1)$, $(r_0,r_1)$ $\in {}_1D_1$. Therefore, $ \partial$ is a derivation and $({}_1D_1, \partial)$ is a differential algebra. It can be exploited to get an automatic way to compute derivatives. Let us consider two functions $f$ and $g$ with their derivatives: their values at the origin, that is the two vectors $(f(0),f’(0))$, $(g(0),g’(0))$, can be considered as two pairs belonging to $_1D_1$. If the derivative of the product $fg$ has to be computed, it is sufficient to look at the second component of the derivative. Defining the operation $[]$ from the space of differential functions to $_1D_1$ via

\[[f] = (f(0),f'(0))\]

we have

\[\begin{gather} [f+g]=[f]+[g],\\ [f\cdot g]=[f]\cdot[g], \\ [1/g]=1/[g]. \end{gather}\]

As a consequence, it is possible to compute the derivatives of many functions by means of the arithmetic rules on $_1D_1$ and by applying the operator $[]$ to the identity function

\[[x]=(x,1).\]

This is equivalent to extract the coefficient of the first order Taylor expansion of the following identity function: \([x]=(x,1)=x+\delta x\).

Consider the following example:

\[f(x) = x^2.\]

Its Taylor expansion centered at $\bar{x}=1$ is equal to

\[T_f(\delta x)=1+2\delta x + \delta x^2.\]

If the function $f$ is evaluated at $(\bar{x},1)=\bar{x}+\delta x$, we get

\[f((\bar{x},1))=(\bar{x}+\delta x)^2=1+2\delta x+\delta x^2.\]

Thus, differential algebra allows an automatic computation of derivatives, that can be easily implemented and exploited in a computer environment. The introduced one can be extended in order to compute derivatives up to an order $n$ of functions in $v$ variables. Moreover, considering that Taylor expansions often give sufficiently accurate approximations of functions, differential algebra provides a tool that can be exploited in several cases. Indeed, by substituting the classical implementation of real algebra with the implementation of a new algebra of Taylor polynomials, any function $f$ can be expanded in its Taylor polynomial up to an arbitrary order, provided that it is sufficiently regular. As an example, let us consider the generic dynamics, given by

\[\dot x = f(x,t),\]

and suppose we know an initial condition $x_0$ with its uncertainty $\delta x_0$: the goal is to perform a propagation up to a time $t_f$. Initialize the initial condition as a vector in $_nD_v$, such that only the first two coefficients are non zero: $[x_0]=x_0+\delta x_0$. If all the operations of the numerical integration scheme are carried out in the framework of differential algebra the state is approximated, at each step $t_i$, as a Taylor expansion centered at $x_0$. By applying Euler’s scheme, at the first step we have

\[[x_1]=[x_0]+f([x_0])\Delta t,\]

where $[x_1]$ may include several non-zeros coefficients corresponding to high order terms in $\delta x_0$. The result of the final step is the $k-th$ order Taylor expansion of the state $x_f$ in $x_0$ at the final time $t_f$: the solution depends on the initial condition in terms of the $k-th$ order polynomial map $T_{x_0}(\delta x_0)$, where $\delta x_0$ is the displacement from the reference initial condition. The accuracy of the result depends on the expansion order $k$. The main advantage of this approach is that it is possible to consider several initial displacements and the new solutions will be obtained by evaluating the polynomial map without any additional numerical integration. Thus, differential algebra can be used as a tool providing a good compromise between computational accuracy and computational burden. For more details, refer to [Berz, 1999] and [Armellin et al., 2010].

# Differential algebra and initial orbit determination

As previously introduced the observations provide the right ascension $\alpha$ and the declination $\delta$ of a celestial body. Three sets $(\alpha_i,\delta_i)$ at three different epochs, $t_i$ with $i=1,2,3$, are sufficient to apply the Gauss method and to determine the body trajectory. However, each observation is affected by an uncertainty: typically, the observation noise can be assumed to be normally distributed. Differential algebra provides a tool to handle this uncertainty. Once the reference solution has been found in terms of position $r_2$ and velocity vector $v_2$ at the epoch $t_2$, differential algebra can be used to obtain its expression as polynomial function of the observations uncertainties:

\[\begin{gather} r_2 = T_{r_2}(\delta \alpha_1,\delta \delta_1,\delta \alpha_2,\delta\delta_2,\delta \alpha_3,\delta \delta_3)\\ v_2 = T_{v_2}(\delta \alpha_1,\delta \delta_1,\delta \alpha_2,\delta\delta_2,\delta \alpha_3,\delta \delta_3)\end{gather}\]

Thus, any probability density function on $(\delta \alpha_i,\delta\delta_i)$ can be mapped on the state $(r_2,v_2)$ using the polynomials $T_{r_2}$ and $T_{v_2}$.

For more details, refer to [Pirovano et al., 2017].