{"id":735,"date":"2018-06-11T02:32:03","date_gmt":"2018-06-11T00:32:03","guid":{"rendered":"https:\/\/aerospaceresearch.net\/?p=735"},"modified":"2018-06-11T12:18:13","modified_gmt":"2018-06-11T10:18:13","slug":"least-squares-fit-to-cartesian-positions-as-a-function-of-time","status":"publish","type":"post","link":"https:\/\/aerospaceresearch.net\/?p=735","title":{"rendered":"[GSoC2018|OrbitDeterminator|Jorge] Week #3-4 \u2013 Least squares fit to cartesian positions as a function of time"},"content":{"rendered":"
\n

Introduction<\/h2>\n
An orbital fit by least-squares is a technique for computing the set of six Keplerian elements which best describe an orbit with respect to a *whole<\/strong>* set of observations.<\/div>\n
<\/div>\n
As the set of Keplerian elements we choose semimajor axis a<\/em>, eccentricity e<\/em>, time of perigee passage \u03c4\u00a0(i.e., how much time has elapsed since the last perigee), as well as the standard Euler angles which describe the orientation of the orbital plane with respecto to the inertial frame: argument of pericenter \u03c9, inclination I<\/em> and longitude of ascending node \u03a9. We denote by X<\/em> the set of these six Keplerian elements; that is, X<\/em> = (a, e, \u03c4, \u03c9, I<\/em>, \u03a9).<\/div>\n
<\/div>\n
A least-squares fit is simply to look the minimum of a scalar (i.e., real-valued) „cost“ function Q(X)<\/em>. The Q<\/em> function is simply a measure of how well our theoretical predictions, for a given set X<\/em> of orbital elements, depart from observations. In an ideal world, we could find a set of orbital elements such that it perfectly matched the observations and thus Q(X)=0<\/em>. Nevertheless, real-world measures carry uncertainties, and we can only hope to find orbital elements X<\/em> which make Q<\/em> attain its minimum possible value with respect to all the observations we have made.<\/div>\n
The values of the cost function Q(X)<\/em> for a least-squared problem are constructed the following way: take *all<\/strong>* the values you observed (call them v_i<\/code>), and substract from each of those the values you were expecting to measure for a given set of orbital elements X<\/em> (say, V_i(X)<\/code> ). Call each of these differences the residuals, xi_i(X) = v_i-V_i(X)<\/code> . Then, square the residuals, and sum all the squared residuals. The value you obtain following this recipe, is the value of Q(X)<\/em>.<\/div>\n
In practice, the least-squares method finds a local minimum of the cost function by performing an iterative procedure, which needs a good initial „guess“ solution X0<\/code>; otherwise the procedure may converge to non-feasible results. Hence, it is necessary that the „good“ initial guess be provided by a preliminary orbit determination method.<\/div>\n

<\/p>\n

<\/div>\n

Implementation<\/h2>\n
We consider the concrete case of having m<\/em> observations of the cartesian coordinates of the position vector (x,y,z)<\/em> referred to the center of attraction (the Earth for artificial satellites, or the Sun for asteroids, etc.) at a given time t<\/em>. That is, for a time t1<\/em> we have (x1,y1,z1)<\/em>; for t2<\/em> we have (x2,y2,z2)<\/em>; etc.<\/div>\n
For this particular case, we assumed the motion to be Keplerian and we wrote the function orbel2xyz<\/code>, which takes as arguments the time t<\/code>, the mass parameter mu<\/code> (mu<\/code> =G \u00d7<\/em>\u00a0m<\/em>) of the center of attraction, the semimajor axis a<\/code>, the eccentricity e<\/code>, the time of pericenter passage tau<\/code>, the argument of pericenter omega<\/code>, the inclination I<\/code> and the longitude of ascending node Omega<\/code>. The output of orbel2xyz<\/code> is the expected value at time t<\/code> of the cartesian positions x<\/code>, y<\/code>, z<\/code>, for a given set of orbital elements. The cartesian position vector is returned as a 1-dim numpy.array<\/code> of length 3<\/code>. Is is constructed by first computing, from the time, the corresponding true anomaly by inverting the Kepler equation n(t-\u03c4)=E-e\u00a0\u00b7 sin(E)<\/em> using a modified version of Newton’s method, as described by Murray and Dermott in their book „Solar System Dynamics“; our implementation is focused on both accuracy and speed and is able to invert Kepler’s equation in 5 iterations to almost machine accuracy, even for high eccentricities (for details, see `time2truean` and related functions). Then, as a second step, we compute the position of the observed object in its orbital frame using standard formulae (this was done in function xyz_orbplane_<\/code>). Finally, using the Euler angles omega<\/code> , I<\/code>\u00a0and Omega<\/code>, we rotate from the orbital frame to the inertial frame using the orbplane2frame_<\/code> function. This allows us to compute our predicted values for the observations, which were denoted by V_i(X)<\/code>\u00a0in the previous section.<\/div>\n
The next step is to retrieve the observed values v_i<\/code> \u00a0from a data file. We do this using the numpy.loadtxt<\/code>\u00a0function.<\/div>\n
The vector of residuals xi_i<\/code> \u00a0is computed via the res_vec<\/code> \u00a0function, which takes as input the set of orbital elements x<\/code>\u00a0as well as a 2-dim numpy.array<\/code> \u00a0called my_data<\/code> , as well as the mass parameter of the Earth in appropriate units, my_mu_Earth<\/code> . Note that if the lengths are, for example, in meters and the times in seconds, then the Earth’s mass parameter has to be given in meters\u00b3\/seconds\u00b2.<\/div>\n
Next, the cost function Q<\/em> is automatically optimized using the scipy.optimize.least_squares<\/code> \u00a0function, taking as input the res_vec<\/code> \u00a0function, an initial guess for the orbital elements x0<\/code> , the complementary arguments for the res_vec<\/code> \u00a0function which in this particular case are the input data data<\/code> \u00a0and the mass parameter of Earth mu_Earth<\/code> . First I thought of using the scipy.optimize.minimize<\/code> \u00a0function; but SciPy’s least_squares<\/code> \u00a0function was chosen because of its access to the Levenberg-Marquardt algorithm, which is quite robust for this kind of problems. The method for nonlinear least squares fitting was the Levenberg-Marquardt algorithm. The output of the least_squares<\/code> \u00a0function is a data structure labeled as Q_ls<\/code>, which among other things, provides the solution to the least-squares problem for the given observational data and associated parameters.<\/div>\n

Results<\/h2>\n
Some results are presented using the determine_kep<\/code> \u00a0method from ellipse_fit<\/code> . For the example file orbit.csv<\/code> , we see an improvement of more than 25% with respect to the initial orbit. Even though the orbital elements are corrected by relatively small quantities, it is sufficient to get a 25% better solution for the orbital elements, in the sense that the cost function attains a lower value with respect to the initial set of elements.<\/div>\n<\/div>\n

<\/p>\n

<\/p>\n

<\/div>\n
\n
\"\"
Least-squares fitting of (t,x,y,z) observations.<\/figcaption><\/figure>\n<\/div>\n
<\/div>\n
\n

Discussion and future work<\/h2>\n
Using well established tools, an implementation of a least-squares orbital determinator for the case of t,x,y,z<\/em> data was implemented. In preliminary tests, we see improvements of around 25% with respect to existing methods. Although execution times are around a couple of seconds, the obtained improvement helps increase the accuracy of the fitted orbital elements, which is essential to orbit propagation.<\/div>\n
Furthermore, although the particular scenario of having observed values of t,x,y,z<\/em> might not be of general use, the tools developed here will help attain least-squares orbit determinators for objects which are observed with optical instruments (telescopes) and range\/range-rate radar astrometry. Indeed, using the code what has been done in this work, it will be feasible to construct an orbit fitter which will handle optical observations. To achieve this, the immediate solutions which must be provided to handle, for example, the optical observations case are the following:<\/div>\n<\/div>\n
<\/div>\n