# Maximum a posteriori estimation through simulated annealing for binary asteroid orbit determination

Irina Kovalenko, Radu Stoica, Nikolay Emelyanov

► To cite this version:

Irina Kovalenko, Radu Stoica, Nikolay Emelyanov. Maximum a posteriori estimation through simulated annealing for binary asteroid orbit determination. Monthly Notices of the Royal Astronomical Society, 2017, 471 (4), pp.4637-4647. 10.1093/mnras/stx1899. hal-02239040

HAL Id: hal-02239040

https://hal.science/hal-02239040v1

Submitted on 19 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Maximum a posteriori estimation through simulated annealing for binary asteroid orbit determination

Irina D. Kovalenko,<sup>1,2★</sup> Radu S. Stoica<sup>1,3</sup> and Nikolay V. Emelyanov<sup>4,1</sup>

Accepted 2017 July 25. Received 2017 July 19; in original form 2016 July 24

## ABSTRACT

In this paper, we consider a new method for the binary asteroid orbit determination problem. The method is based on the Bayesian approach with a global optimization algorithm. The orbital parameters to be determined are modelled through an a posteriori distribution made of a priori and likelihood terms. The first term constrains the parameter space and it allows the introduction of available knowledge about the orbit. The second term is based on given observations and it allows us to use and compare different observational error models. Once the a posteriori model is built, the estimator of the orbital parameters is computed using a global optimization procedure: the simulated annealing algorithm. The maximum a posteriori (MAP) techniques are verified using simulated and real data. The obtained results validate the proposed method. The new approach guarantees independence of the initial parameter estimation and theoretical convergence towards the global optimization solution. It is particularly useful in the following situations: whenever a good initial orbit estimation is difficult to obtain, whenever observations are not well sampled and whenever the statistical behaviour of the observational errors cannot be said to be Gaussian-like.

**Key words:** methods: statistical – celestial mechanics – minor planets, asteroids: general.

## 1 INTRODUCTION

The estimation of orbital parameters from observations is a classical inverse problem of celestial mechanics and it has been developed extensively for single asteroids. For binary asteroid mutual orbits, the problem is made complicated by the larger number of unknown parameters (as the mass of the system is usually unknown) and also by insufficient observational data.

Several methods have been developed for the determination of binary objects' mutual orbits. The classical least-squares (LS) orbit fitting to observations has been used for satellites of asteroids (Descamps 2005) as well as for satellites of planets (Emelyanov 2005). For trans-Neptunian binaries (TNBs), Noll et al. (2004) used the downhill simplex Nelder-Mead algorithm (Nelder & Mead 1965). Vachier, Berthier & Marchis (2012) developed a statistical method for mutual orbit fitting using a genetic-based algorithm. All these methods are extremely efficient from a numerical point of view for the orbit determination problem, and they give accurate results. However, they require good initial estimations of

the involved model parameters, otherwise the algorithms converge towards local optima of the fitting parameters.

Bayesian a posteriori probability density formulation into orbit determination was introduced by Muinonen & Bowell (1993). Their approach was implemented in the so-called statistical ranging method by Virtanen, Muinonen & Bowell (2001), which uses Monte Carlo techniques for asteroid heliocentric orbit determination. The statistical ranging method constrains a volume of orbits in the orbital-element phase space, using the sampling in the observational space. Namely, each sample orbit is produced from two randomly chosen observations, introducing random deviations in the observations. Finally, a large set of orbits compatible with the observations is produced, while the six-dimensional element space is constrained. Oszkiewicz et al. (2009) developed the orbit ranging method that uses a Markov chain Monte Carlo (MCMC) method based on the Metropolis-Hastings (MH) algorithm. Another MCMC method, called the virtual-observation MCMC method, was proposed by Muinonen et al. (2012). In this method, the sampling is done in the parameter space instead of the orbital space. Thus, each sample orbit is generated from virtual observations with the downhill simplex Nelder-Mead algorithm.

For binary mutual orbit determination, Grundy et al. (2008) used the aforementioned Monte Carlo techniques (Muinonen &

<footnote n="1"><sup>1</sup>IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Lille, 77 av. Denfert-Rochereau, F-75014, Paris, France</footnote>

<footnote n="2"><sup>2</sup>Space Research Institute, Russian Academy of Sciences, 13 Profsoyuznaya ul. 84/32, 117342 Moscow, Russia</footnote>

<footnote n="3"><sup>3</sup>Université de Lorraine, Institut Elie Cartan de Lorraine, F-54506 Vandoeuvre-lès-Nancy Cedex, France</footnote>

<footnote n="4"><sup>4</sup>M. V. Lomonosov Moscow State University – Sternberg Astronomical Institute, 13 Universitetskij prospect, 119992 Moscow, Russia</footnote>

<footnote>* E-mail: irina.kovalenko@obspm.fr</footnote>
Bowell 1993; Virtanen et al. 2001) in addition to the downhill simplex algorithm (Nelder & Mead 1965). Hestroffer, Vachier & Balat (2005) and Oszkiewicz, Hestroffer & David (2013) developed techniques for binary asteroids based on Monte Carlo sampling with the Thiele-Innes method (Aitken 1964). The Thiele-Innes method, initially proposed for the determination of the orbits of binary stars, requires the orbital period and three observations in the same tangent plane. Oszkiewicz et al. (2013) use these parameters as the sampling parameters in the MH algorithm, so that the orbit sample is produced by varying three observations in the same tangent plane, randomly selected from the whole set of observations, and orbital period. Even though the approach of Oszkiewicz et al. (2013) has already been successfully applied to a few binaries, generally the use of the Thiele-Innes method might be inefficient for binary asteroids. The inherent parallax, because of the relative motion between the observer and the target, cannot be neglected, even for observations on close tangent planes, and the random selection of three observations in the same tangent plane is not usually possible for a given data set.

The method we propose in this paper continues and develops the previous ideas of Bayesian modelling. The a posteriori model contains two terms: the likelihood and the a priori. The likelihood term relies on the orbital parameters with the observations. In this work, the likelihood is built so that different observational error distributions can be considered, whereas most methods usually assume only Gaussian observational errors. The a priori allows us to introduce available knowledge about the orbit and to constrain the parameter space. By analogy with the optimization problem, it plays the role of a regularization term. This approach is useful in the case of newly discovered binary asteroids for preliminary orbit determination, when the initial parameter estimation is difficult. Once the likelihood and prior are built, the estimate of orbital parameters is given by the MAP method. The MAP estimate is computed using the simulated annealing (SA) algorithm. This algorithm is a global optimization method that guarantees convergence towards the global optimum of a considered function, independently of the initial state. The orbit sampling in the SA method is performed by the MH algorithm through the variation of orbital parameters. The advantage of such an approach, over the sampling through observations, is the avoidance of parameter computation from observations (e.g. by the Thiele-Innes method) for each orbit sample.

The structure of the paper is as follows. In Section 2, we describe the dynamical model with the observations of binary asteroids. Then, we present the problem in the form of a non-linear regression model in Section 3. We describe the LS estimator and the new estimator that we propose, associated with the SA algorithm, in Sections 4 and 5, respectively. We describe the validation of the new method on simulated and real observations in Section 6. Finally, we give conclusions and perspectives in Section 7.

## 2 DYNAMICAL MODEL

### Binary system

A binary asteroid is a system of two asteroids orbiting their common centre of mass. Here, it is assumed that the motion of such systems follows Kepler's laws, and the attractive force in the binary system is a gravitational force expressed by Newton's law.

### Observations and coordinates

Unless the relative masses of two components are known, it is impossible to determine the centre of mass and therefore the orbits around the centre. However, we are able to derive the Keplerian relative orbit of one component described about the other, which is assumed to remain stationary at the focus.

The motion that we observe is not the true motion but its projections on to planes perpendicular to the line of sight. The observation gives us the apparent position of one asteroid, called the companion or secondary, with respect to another, called the primary, which is assumed to be at the focus. Every complete observation of a visual binary asteroid supplies the time of observation *t* and two coordinates on the sky plane. Let us set rectangular coordinates $x_t$, $y_t$ measured on the tangent plane. The position of the tangent plane is determined by $\alpha$ and $\delta$, right ascension and declination, respectively, at a given time *t*. Measured in arcseconds, the tangent coordinates are related to $\alpha$, $\delta$ through the equation $x_t = (\alpha_2 - \alpha_1)\cos\delta_1$, $y_t = \delta_2 - \delta_1$, where $(\alpha_1, \delta_1)$ and $(\alpha_2, \delta_2)$ are related to the primary and secondary, respectively.

We introduce a coordinate frame, related to the tangent plane and centred at the primary asteroid. The *y*-axis is directed to the north, the *x*-axis to the east, and the *z*-axis is normal to the tangent plane and directed away from the observer. In this frame, the relative position of the secondary asteroid with respect to the primary one is related to observed positions through $x = Rx_t$ and $y = Ry_t$, where *R* is the distance from the observer to the asteroid and the coordinates $x_t$ and $y_t$ are transformed from arcsec to radians.

Let us introduce an equatorial coordinate system ($x_E$, $y_E$, $z_E$) referenced to the epoch J2000 with the centre at the primary asteroid. Thus, the rotation matrix from the tangent frame (*x*, *y*, *z*) to the equatorial frame ($x_E$, $y_E$, $z_E$) is

$$S(\alpha, \delta) = \begin{pmatrix} -\sin\alpha & -\cos\alpha \sin\delta & \cos\alpha \cos\delta \\ \cos\alpha & -\sin\alpha \sin\delta & \sin\alpha \cos\delta \\ 0 & \cos\delta & \sin\delta \end{pmatrix}$$

The topocentric distance *R* as well as the right ascension $\alpha$ and declination $\delta$ of an asteroid at the moment of observation can be found from the ephemerides of asteroid heliocentric motion. They can be treated here as known values.

### Time delay of the light

The time delay of the light between the observation and asteroid position has to be taken into account and calculated by the linear equation $t = t_0 - R/c$, where $t_0$ is the time of observation, *R* is the distance between the observer and target and *c* is the light velocity.

### Orbital parameters

We consider the Keplerian orbit of the secondary asteroid, with the primary being located in the focus of the ellipse. This orbit is described by seven parameters (*a*, *e*, *i*, $\Omega$, $\omega$, $\tau$, *P*): the six Keplerian elements – semimajor axis *a*, eccentricity *e*, inclination *i*, longitude of the ascending node $\Omega$, argument of pericentre $\omega$, time of pericentre passage $\tau$ – and the orbital period *P*, which is also an independent parameter under the condition that the mass of the system is also unknown. Otherwise, it can be derived using the third Kepler's law. The angles *i*, $\Omega$ and $\omega$ are referenced to the J2000 equatorial frame.
## Coordinates from orbital parameters

The direct problem is to calculate the position of an asteroid at a given time from a known orbit. When the elements $(a, e, i, \Omega, \omega, \tau, P)$ are known, the apparent position coordinates $(x, y)$ at a given time $t$ are derived as follows. The mean anomaly:

$$M = n(t - \tau),$$

where $n = 2\pi/P$. The eccentricity $E$ is determined from the Kepler equation: $M = E - e\sin E$.

The spherical coordinates – the radius vector $r$ and the true anomaly $\nu$ – are calculated using the following equations:

$$r = a(1 - e \cos E); \quad \tan \nu = \frac{\sqrt{1 - e^2} \sin E}{\cos E - e}.$$

The coordinates in the equatorial frame are

$$x_E = r(\cos u \cos \Omega - \sin u \sin \Omega \cos i),$$

$$y_E = r(\cos u \sin \Omega + \sin u \cos \Omega \cos i),$$

$$z_E = r \sin u \sin i,$$

where $u = \nu + \omega$ is argument of latitude.

The final $(x, y)$ apparent position coordinates are

$$x = -x_E \sin \alpha + y_E \cos \alpha,$$

$$y = -x_E \cos \alpha \sin \delta - y_E \sin \alpha \sin \delta + z_E \cos \delta.$$

Once the $(x, y)$ coordinates are calculated for a given set of orbital parameters, the observed minus calculated residuals (O – C) can be computed for each observation.

## 3 REGRESSION MODEL

The regression model construction needs to relate observations and theoretical positions to each other. Let each observation at a given time $t^{(k)}$ consist of two measured values $x_k^o, y_k^o$, $k = 1, 2, \dots, N$, where $N$ is the number of observations. Let us denote

$$\varphi_{2k-1} = x_k^o,$$

$$\varphi_{2k} = y_k^o,$$

$$t_{2k-1} = t_{2k} = t^{(k)}.$$

Thus, the set of observations can be described by the following vector:

$$\varphi = (x_1^o, y_1^o, x_2^o, y_2^o, \dots, x_N^o, y_N^o)^T \\ = (\varphi_1, \varphi_2, \varphi_3, \dots, \varphi_{2N})^T.$$

With a given theoretical model and time $t$ of observation, the measured values $x$ and $y$ can be computed. The theory (see Section 2) provides two functions:

$$x = \psi_x(t, \theta), \quad y = \psi_y(t, \theta),$$

with $\theta = \{\theta_j\}, j = 1, \dots, 7$ model parameters.

For a set of moments of time $t = \{t_i\}$ ($i = 1, 2, \dots, 2N$), where $t_i = t^{(k)}, k = i/2$ when $i$ is even and $k = (i + 1)/2$ when $i$ is odd, we denote

$$\Psi = (\psi_1, \psi_2, \psi_3, \dots, \psi_{2N})^T,$$

where

$$\psi_i = \begin{cases} \psi_x(t_i, \theta), & \text{when } i \text{ is odd;} \\ \psi_y(t_i, \theta), & \text{when } i \text{ is even.} \end{cases}$$

The so-called observational equation, which relates observations and theoretical positions to each other, can be expressed as

$$\varphi = \Psi(\theta) + \varepsilon, \quad (1)$$

where $\theta$ is a vector of seven unknown parameters that describe the orbit (see Section 2) and $\varepsilon$ describes observational and theoretical errors, which are also unknown, $\varepsilon = (\varepsilon_{x1}, \varepsilon_{y1}, \dots, \varepsilon_{xN}, \varepsilon_{yN})^T$. Equation (1) represents a non-linear regression model. Statistical inference about unknown parameters $\theta$ and $\varepsilon$ can be drawn from $\varphi$. Two methods for the problem of finding $\theta$ are described here: the classical LS estimator and the new MAP estimate based on the SA algorithm.

## 4 LEAST-SQUARES METHOD

A commonly used method in celestial mechanics is to fit orbital parameters using the algorithm of the LS approach. Here, the application of this method is described for the determination of the orbital parameters of binary asteroids.

Let us use the same designations as in the observational equation (1). If observational and theoretical errors are neglected, a conventional observational equation can be written as

$$\varphi = \Psi(\theta) \quad \text{or} \quad \varphi_i = \psi_i(\theta) \quad (i = 1, 2, \dots, 2N),$$

where $\theta = (a, e, i, \Omega, \omega, \tau, P) = \{\theta_j\}, j = 1, \dots, 7$ is the vector of the true parameter values. Let $\theta^0 = \{\theta_j\}, j = 1, \dots, 7$ be the vector of initial parameters. Thus, the correction $\Delta\theta_j$ between each initial and true $j$th parameter is $\Delta\theta_j = \theta_j - \theta_j^0$, and the conventional observational equation is

$$\varphi_i = \psi_i(\theta^0 + \Delta\theta), \quad (i = 1, 2, \dots, 2N), \quad (2)$$

where $\Delta\theta = \{\Delta\theta_j\}, j = 1, \dots, 7$. The LS method assumes that the estimation of the preliminary orbital parameters is close enough to the final solution that it can be used to adjust parameters rather than to determine them. Hence, this allows us to assume that $\Delta\theta$ is small and to expand equation (2) to the Taylor series. Let us retain only first-order terms of $\Delta\theta_j$, so

$$\Delta\varphi_i = \sum_{j=1}^{7} J_{i,j} \Delta\theta_j \quad (i = 1, 2, \dots, 2N), \quad (3)$$

where we denote $\Delta\varphi_i = \varphi_i - \varphi_i^{c(0)}, \varphi_i^{c(0)} = \psi_i(\theta^0)$ and $J_{i,j} = \partial\psi_i(\theta^0)/\partial\theta_j$. The approximation (3) is the so-called conditional equation. These conditional equations are approximate, because the errors of the theoretical model as well as the observational errors are neglected here, as are all terms higher than first-order terms in the Taylor series. Once the adjustments $\Delta\theta_j$ are derived, they are added to the initial parameters. This parameter determination method is called differential adjustment. The adjustments can be repeated many times, and, if the process converges (i.e. if adjustments become smaller), then the algorithm can be stopped as soon as the adjustments are smaller than the assumed uncertainties.

At each step of parameter adjustment, the residuals of the conditional equation (3) are given by

$$s_i = \Delta\varphi_i - \sum_{j=1}^{7} J_{i,j} \Delta\theta_j \quad (i = 1, 2, \dots, 2N).$$

There are many algorithms of parameter fitting. Here we use the least-squares fitting that minimizes the sum of the squared residuals. To minimize the sum of squares of $s_i$, the gradient equation is set
to zero and is solved for $\Delta\theta_j$. We obtain seven simultaneous linear equations, which are called normal equations and are written in matrix notation as

$$ (J^T J) \Delta\theta = J^T \Delta\varphi. $$

## 4.1 Weighted least-squares method

The classical LS method can be generalized by the weighted least-squares (WLS) method. This method considers the case when observations have different variances by introducing weights. Thus, if the measurements are uncorrelated, the normal equations should be modified as

$$ (J^T W J) \Delta\theta = J^T W \Delta\varphi, $$

where $W$ is a diagonal weighted matrix, with $W_{i,i} = 1/\sigma_i^2$, where $\sigma_i^2$ is a variance estimate of the $i$th measure.

# 5 BAYESIAN MODELLING AND SIMULATED ANNEALING

According to the Bayesian framework, a posteriori probability density function of the unknown parameters $\theta$ given the observed data $\varphi$ can be written as

$$ p(\theta|\varphi) \propto \mathcal{L}(\theta, \varphi) p(\theta), \quad (4) $$

where $\mathcal{L}(\theta, \varphi)$ and $p(\theta)$ are the likelihood and a priori terms, respectively. The role of the likelihood term is to find the ‘best’ orbital parameters to fit the observations. Nevertheless, using the likelihood term only will give an ill-conditioned problem. This drawback is solved by using a regularization term – the prior density.

Hence, under these considerations, the parameter estimator we propose is the MAP, given by

$$ \hat{\theta} = \arg\max_{\theta \in \Theta} p(\theta|\varphi), \quad (5) $$

where $\Theta \subset \mathbb{R}^7$ is the parameter space.

## 5.1 Likelihood model: conditional data term

Hereafter, for the sake of simplicity and in agreement with the observational equation (1), it is convenient to consider $N$-dimensional vectors instead of $2N$. Thus, $\varphi = \{\varphi_i\}, (i = 1, 2, \dots, N)$, where $\varphi_i = (x_i^O, y_i^O)$ is the observed position of the secondary asteroid with respect to the primary, $t = (t_1, t_2, \dots, t_N)$, $t_i$ is the time of observation, $\varepsilon = (\varepsilon_{x1}, \varepsilon_{y1}; \varepsilon_{x2}, \varepsilon_{y2}; \dots; \varepsilon_{xN}, \varepsilon_{yN})^T$ are the errors and $\psi = \{x_i, y_i\}, (i = 1, 2, \dots, N)$ is the set of theoretical positions. As before, the vector of unknown parameters is $\theta = (a, e, i, \Omega, \omega, \tau, P)$.

For the construction of our model, no particular assumption is favoured. Nevertheless, it is reasonable to assume that the observations are fitted as closely as possible by the optimal orbital parameters. One way to obtain this fitting is to allow the likelihood to depend on the distances from observed points to their associated calculated points. The calculated point is the point position computed using the current orbital parameters at the same time $t$ as the given observed point.

The likelihood function in equation (4) is defined by

$$ \mathcal{L}(\theta, \varphi) = \prod_{i=1}^{N} p(\varphi_i|\theta), $$

where we set

$$ \begin{aligned} p(\varphi_i|\theta) &= \exp \left[ - \left( |x_i^O - x_i^C|^l + |y_i^O - y_i^C|^l \right)^{k/l} \right] \\ &= \exp[-\varepsilon_i(l)^k]. \end{aligned} $$

Here, $(x_i^O, y_i^O)$ and $(x_i^C, y_i^C)$ are the coordinates of observations and the computed positions at the same associated times, respectively. The functions $\varepsilon_i(l)$ can be interpreted as the equivalent of the residuals in a regression model. For instance, whenever $(k, l) = (1, 1)$, the residuals tend to exhibit a Laplacian character, while they show a Gaussian character for $(k, l) = (2, 2)$.

Furthermore, whenever the uncertainties of observations are available, the likelihood is written as

$$ \mathcal{L}(\theta, \varphi) = \exp \left[ - \sum_{i=1}^{N} \varepsilon_i(l)^k \right], \quad (6) $$

while the residuals are given by

$$ \varepsilon_i(p) = \left[ \left( |x_i^O - x_i^C| / \sigma_{x,i} \right)^l + \left( |y_i^O - y_i^C| / \sigma_{y,i} \right)^l \right]^{1/k/l}, \quad (i = 1, \dots, N), $$

where $\sigma_{x,i}$ and $\sigma_{y,i}$ are the estimated uncertainties for a given observation.

In this work, the following four likelihood models (6) are compared:

model 1, $(k, l) = (1, 1)$, $\sigma_{x,i} = \sigma_{y,i} = 1, (i = 1, 2, \dots, N)$;

model 2, $(k, l) = (2, 2)$, $\sigma_{x,i} = \sigma_{y,i} = 1, (i = 1, 2, \dots, N)$;

model 3, $(k, l) = (1, 1)$, $\sigma_{x,i}, \sigma_{y,i} (i = 1, 2, \dots, N)$ are the given uncertainties, estimated for each observation;

model 4, $(k, l) = (2, 2)$, $\sigma_{x,i}, \sigma_{y,i} (i = 1, 2, \dots, N)$ are the given uncertainties, estimated for each observation.

## 5.2 A priori model: regularization term

The a priori model allows the introduction of available knowledge regarding the joint probability distribution of the orbital parameters. For instance, if correlations between the different parameters are previously known, they can be introduced through the prior model. This term also defines the restriction of the parameter space to those of the considered celestial bodies.

Whenever no particular information is available, a prior should have minimal influence on the inference. A possible choice is to use the Jeffreys non-informative priors (Jeffreys 1946). Jeffreys noted that for independent parameters the prior should be treated separately. Thus, in our problem, it is convenient to use

$$ p(\theta) = \prod_{i=1}^{7} p(\theta_i), $$

where $\theta_i$ is one of the seven parameters describing the orbit.

Following the Jeffreys principle of non-informative prior choice when the parameter space is a bounded interval (Jeffreys 1973), we set the independent prior distributions for each parameter to be uniform. These distributions were defined over the largest interval of possible values corresponding to each parameter.

### Semimajor axis

The lower limit on the semimajor axis can be derived from observations. Let $\rho$ be the observed distance between two components in a binary system. As a projection on to the tangent plane, the maximal
value $\rho$ is never greater than the sum of the semimajor axis $a$ and the distance $c$ from the focus to the centre of the orbit ellipse: $\rho_{max} \le a + c$, where $c = ae$. In the case of a very elongated orbit, the eccentricity is $e \approx 1$, and thus $\rho_{max} \le 2a$. This allows us to set the lower limit on the semimajor axis: $a_{min} = 0.5\rho_{max}$. While the lower limit can be derived, an approximation is needed for the upper limit. Namely, the a priori of the semimajor axis is set to be

$$a \sim \mathcal{U}[a_{min}, a_{max}],$$

where $a_{max}$ is the maximum semimajor axis among known binaries in the same dynamical group (e.g. near-Earth, main-belt, trans-Neptunian object).

### Eccentricity, longitude of the ascending node and argument of periapsis

In the absence of a priori information about the eccentricity $e$, the longitude of the ascending node $\Omega$ and the argument of periapsis $\omega$, we assume that they follow uniform distributions in the interval of all possible values:

$$e \sim \mathcal{U}[0, 1], \quad \Omega \sim \mathcal{U}[0, 2\pi], \quad \omega \sim \mathcal{U}[0, 2\pi].$$

### Inclination

When the topocentric direction of the asteroid is constant, two mirror orbits are possible: direct and retrograde. The first has an inclination between 0 and 90 degrees, and the second between 90 and 180 degrees. Thus, the problem can be split for finding two solutions that correspond to the following a priori:

$$i \sim \mathcal{U}[0, \frac{\pi}{2}], \quad \text{direct motion}$$

or

$$i \sim \mathcal{U}[\frac{\pi}{2}, \pi], \quad \text{retrograde motion.}$$

In the case of distant asteroids or TNBs, the topocentric direction can be almost constant. In general case when the orbit is observed from different sides (with different topocentric directions) the inclination value can be distinguished from observations.

### Orbital period

The interval of possible values for the orbital period cannot be constrained analytically from observations and/or the dynamical model. In order to set a prior for this parameter, approximations were made based on the properties of the asteroid dynamical group. We denote the a priori distribution

$$P \sim \mathcal{U}[P_{min}, P_{max}],$$

where $P_{min}$ and $P_{max}$ are the minimal and maximal values of the orbital period among known binaries in the same dynamical group. Additionally, the average asteroid, $a_{min}$ and the third Kepler's law can be used in order to restrict the prior interval.

### Time of periapsis passage

The time of periapsis passage $\tau$ can be bounded by the moment of the first observation $t_{min}$. The a priori distribution for $\tau$ is uniform in the range equal to the period

$$\tau \sim \mathcal{U}[t_{min}, t_{min} + P_{max}],$$

where $P_{max}$ is the upper orbital period a priori limit.

## 5.3 Optimization algorithm: the simulated annealing

The SA algorithm has been built to compute the MAP estimate (5). This algorithm is a general global optimization technique that works by iteratively simulating $p(\varphi|\theta)^{1/T}$ while slowly cooling the temperature $T$ (Kirkpatrick, Gelatt & Vecchi 1983). The principle of the method is the following. At high temperatures, all the configurations are accepted, and the simulated probability distribution is equivalent to the uniform distribution over the entire configuration space. While the temperature decreases, the algorithm chooses those states that tend to maximize the considered probability distribution. If the temperature decreases slowly enough, then the algorithm makes its choices while getting out from the local optima. At convergence, when temperatures $T \to 0$, the algorithm is frozen at the desired global maximum. In theory, the algorithm converges weakly to the uniform distribution over the subspace of configurations that maximize the probability density of interest (Geman & Geman 1984; Stoica, Gregori & Mateu 2005).

Two ingredients are needed to set up such an algorithm. The first is a sampling algorithm from the considered probability distribution, and the second is a cooling schedule for the temperature parameter. For the sampling algorithm, we chose to build MH dynamics (Hastings 1970). This algorithm converges theoretically after an infinite number of iterations, and the speed of its convergence depends on the proposal densities used. The initial set of parameters $\theta^{(0)}$ can be chosen according to the a priori distribution. Being in the state $\theta^{(i-1)}$, a new parameter value $\theta'$ is proposed using a uniform proposal density for each parameter

$$Q(\theta'|\theta^{(i-1)}) = \prod_{j=1}^{7} \frac{1\{\theta'_{j} \in [\theta_{j}^{(i-1)} - \Delta\theta_{j}, \theta_{j}^{(i-1)} + \Delta\theta_{j}]\}}{2\Delta\theta_{j}}. \quad (7)$$

On the one hand, if $\Delta\theta_j$ is too high, the new proposed state risks being always rejected. On the other hand, if $\Delta\theta_j$ is too low, the algorithm does not travel fast enough through the configuration space. For our problem, the values for $\Delta\theta_j$ were found after several trials and errors. Here, because the parameter space is compact in $\mathbb{R}^d$, these choices for the proposal densities guarantee that the algorithm converges with equal speed, independently of initial state (Roberts & Smith 1994; Roberts & Tweedie 1996).

Geman & Geman (1984) and Stoica et al. (2005) proved the existence of logarithmic cooling schedules that guarantee the convergence of their proposed SA algorithms associated with their corresponding simulation dynamics. Establishing such a result for any probabilistic model and simulation dynamics is a non-trivial task. From a more practical perspective, a common choice for the cooling schedule is

$$T_i = T_{max} c^i, \quad (8)$$

where $T_{max}$ is the initial temperature and $c$ is a constant in the interval $[0.95, 1)$. As an option, it might be preferable to slow down this exponential cooling schedule, while reducing the correlation between samples. This can be done by reducing the temperature every $m$ iterations, where $m$ is the so-called decorrelation step.

## 5.4 Summary of the algorithm

The implemented SA algorithm consists of the following steps.

### Step 1. Initialization.

(a) Choose an initial value for the orbital parameters $\theta_0 \in \Theta$. The parameters are randomly chosen according to the prior distributions.
(b) Set an initial temperature $T_{max}$.

(c) Set a counter $i = 1$.

(d) Set stop conditions: a limiting number of successive temperatures $n_{max}$ when no new proposal orbit is accepted (the system is called 'frozen'; Li 2009) and a limiting number of iterations $i_{max}$.

(e) Set a decorrelation step $m$.

Step 2. For each iteration $i$, do the following.

(a) Generate a candidate set of orbital parameters from the proposal distribution

$$ \theta' \sim Q(\theta'|\theta^{(i-1)}), $$

where $Q(\theta'|\theta^{(i-1)})$ is given by equation (7).

(b) Compute the acceptance probability:

$$ \alpha = \min \left\{ 1, \left[ \frac{p(\theta'|\varphi)}{p(\theta^{(i-1)}|\varphi)} \right]^{1/T_{i-1}} \right\}. $$

(c) The new orbit is accepted with probability $\alpha$.

In practice, a random number $\alpha_r$ is generated in the range [0, 1] and is compared with $\alpha$. If $\alpha_r < \alpha$, the new orbit is accepted $\theta^{(i)} = \theta'$; otherwise, it is rejected and the last accepted orbit is set as the current one $\theta^{(i)} = \theta^{(i-1)}$.

(d) If $(i \mod m) = 0$, set a new temperature $T_i = h(i, T_{max})$, where $h$ is the cooling schedule (8).

(e) Set $i = i + 1$.

Step 3. The algorithm stops if no new proposed orbit is accepted during $n_{max}$ successive changes of the temperature, or if it achieves $i = i_{max}$.

## 5.5 Uncertainty estimation

In order to obtain a statistical evaluation of the result, the algorithm 5.4 was performed 100 times using different initial sets of orbital parameters $\theta_0$, which were chosen following the prior distribution. The initial temperature $T_{max}$ was the same for all of these situations. Thus, we obtain 100 orbits from which we retain the one (described by the set of orbital parameters $\hat{\theta}$) that produces the greatest a posteriori value.

The uncertainty estimation of the parameters is given by the empirical variance and an inter-quantile interval. The variance for each parameter is estimated from the obtained 100 orbit samples. Then the uncertainty is given by two estimated standard deviations $\pm 2\sigma_{\theta}$. The length of the chosen inter-quantile interval of level 95 per cent for each parameter is computed using the difference of the corresponding empirical quantiles $q(0.975)$ and $q(0.025)$. In addition, the length of this interval can be compared with $4\sigma_{\theta}$ as a verification of the central limit theorem.

The approach of uncertainty estimation through inter-quantile intervals is more robust to those values that tend to be either too small or too large. This is not the case for the variance estimate, as this estimate depends on the sum of all of the considered values. Under these circumstances, for the ephemeris uncertainty estimation we use only the inter-quantile interval. The uncertainty for a predicted position at a given time $t$ is estimated by a 95 per cent inter-quantile interval as follows. Once the sample of orbits is obtained, we calculate the sample of positions at time $t$. For the obtained sample, the 95 per cent inter-quantile interval is obtained using the empirical quantiles $q(0.025)$ and $q(0.975)$ for the position coordinates. If the coordinates of the observed position at the corresponding time $t$ are not inside the inter-quantile interval, the hypothesis that the

observed position might be a realization of the estimated model is rejected, with a p-value lower than 5 per cent.

# 6 PRACTICAL ORBIT COMPUTATION

## 6.1 Implementation

The input data consist of a set of observations, the likelihood and the a priori models, and the SA configurations. For each observation at time $t_i$ ($i = 1, 2, ..., N$) the following data should be provided: the coordinates $x_i, y_i$ of the relative position of the secondary with respect to the primary, the topocentric distance $R_i$, and the right ascension $\alpha_i$ and declination $\delta_i$ of the primary in the geocentric frame of the J2000 epoch.

The SA algorithm parameters (the proposal distributions and the temperature schedule) were set after several trials and errors. The proposal distribution was chosen to be the uniform distribution with support $\Delta\theta_j = 0.1\Delta\theta_{(prior)j}$, where $\Delta\theta_{(prior)j}$ defines the associate uniform prior. There is a lot of freedom to choose the proposal distribution, guaranteeing the theoretical convergence of the SA, while preserving the uniform ergodicity of the simulated Markov chain at the basis of the optimization algorithm. Nevertheless, choosing the optimal proposal distributions in order to obtain the highest speed of convergence is an open mathematical problem. Therefore, the previously mentioned choice was preferred not only in agreement with the theoretical requirements, but also because of its simplicity. The cooling is performed according to the formula (8), with decorrelation step $m = 50$, coefficient $c = 0.999$ and $T_{max} = 10^7$. The algorithm stops if no proposal orbit is accepted at 100 successive temperatures.

Experiments on simulated and real data sets were conducted in order to verify the proposed method. The considered data sets are related to TNBs and are presented in the Appendix.

In the following, the four likelihood models described in Section 5.1 are analysed. The prior distribution is described in Section 5.2. In the particular case of TNBs, the maximum semimajor axis for the prior is given $a_{max} = 102000$ km (Johnston 2014). The orbital period according to Johnston (2014) ranges from $\approx 0.5$ to 6300 d. Because this interval is large, we restrict the prior and choose $P \sim U[0.5, 1000]$. This estimation is in accordance with the average mass of known TNBs $m \approx 4.5 \times 10^{18}$ kg and the calculated $a_{min}$ (according to the third Kepler's law).

## 6.2 Simulated observations

Ten artificial observations were generated using the model of the binary asteroid motion, described in Section 2, and the following orbital parameters: $a = 10000$ km, $e = 0.5$, $i = 135^\circ$, $\Omega = \omega = 45^\circ$, $\tau = 0$, $P = 30$ d. The right ascension $\alpha$, declination $\delta$ and geocentric distance correspond to those of the TNB Altjira. Here we assume that the difference between the geocentric and topocentric $R$ distances is neglected. The values of coordinates $\alpha, \delta$ and $R$ are set from the ephemerides of its heliocentric motion on the range between JD 245 4000.5 (for the first observation) and JD 245 4057.5 (for the last observation); see Table A1. The simulated observations cover one orbital revolution of the secondary around the primary. The choice of the distant trans-Neptunian object allows us to avoid a significant tangent plane changing during the considered interval of time, and to graphically demonstrate the distribution of observations (see Fig. 1). The 10 generated observations of the relative positions, listed in Table A1, are non-homogeneously distributed
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Simulated observations (black points): $\Delta x$ and $\Delta y$ correspond to relative positions of the secondary with respect to the primary. Lines show a search for the optimal solution during the SA algorithm (only each 100th change of calculated position, from the last 7000 accepted orbits, is shown).</figcaption></figure>

along the apparent orbit (see Fig. 1). We randomly set the observational uncertainties in the range from 0.001 to 0.01 arcsec following the uniform distribution (see Table A1).

## Results

An example of the SA behaviour is shown in Fig. 1. This plot shows how the calculated positions evolve towards the matching of the observations during the MAP search.

For the statistical inference, we repeat the SA algorithm 100 times, giving different random initial sets of orbital parameters (see Section 5.5). The resulting distribution of orbit samples, obtained with the first model of the likelihood, is shown in Figs 2 and 3. The distributions obtained with the other three models are similar to those of the first model and are not graphically shown here. According to the obtained distributions, we retain the mean value, the standard deviation and the 2.5–97.5 per cent quantiles for each orbital parameter. This statistical summary for the four likelihood models is listed in Table 1. The mean values of the obtained distributions (see Table 1) are very close to the original orbital parameters. This result is obtained despite the fact that the search interval for the parameters, given through the a priori term, was chosen to be as large as possible. This also demonstrates how the MAP converges towards the optimal solution.

The uncertainties from models 1 and 3 are greater than those for models 2 and 4, respectively, following both interval estimations – 2σ and the inter-quantile. Still, for all the models, there is an agreement between the length of the inter-quantile interval and $4\sigma_{\theta}$. This indicator, together with the symmetric bell-shape of the parameter distributions (as presented by histograms in Figs 2 and 3), suggests no evidence against un-biasedness and Gaussian statistical behaviour of the obtained results. Because the chosen a priori was uniform for all the parameters, this tends to be consistent with the asymptotical properties of the maximum likelihood estimates.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Resulting distributions of semimajor axis *a*, eccentricity *e*, inclination *i* and longitude of the ascending node $\Omega$ (referenced to the J2000 equatorial frame) obtained for simulated observations with the likelihood model 1. The solid line represents the true parameter value. The dotted lines represent the 2.5 and 97.5 per cent quantiles of the resulting sample.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Resulting distributions of argument of pericentre $\omega$, mean anomaly *M* (at 15th simulation day, referenced to the J2000 equatorial frame), orbital period *P* and the average O – C residual obtained for simulated observations with the likelihood model 1. The solid line represents the true parameter value. The dotted lines represent the 2.5 and 97.5 per cent quantiles of the resulting sample.</figcaption></figure>

In addition, the average sky plane residual for each orbit sample is done using the following formula:

$$ \overline{O-C} = \frac{1}{N} \sum_{i=1}^{N} \left( |x_i^{\text{o}} - x_i^{\text{c}}|^2 + |y_i^{\text{o}} - y_i^{\text{c}}|^2 \right)^{1/2}, \quad (9) $$

The ranges of the average O – C values for the obtained sample of orbits, associated with different likelihood models in the MAP algorithm, are listed in Table 2. The histogram of the sky plane residuals, shown in Fig. 3, indicates the prevalence of small values, and hence a good fit of the model.

To summarize, no remarkable advantage of any likelihood model over others is noticed. The method converges to the optimal orbit equally with the four models. The outcome of different uncertainties does not favour one model over the others.
<table><caption>**Table 1.** Statistical summary of parameter distribution, obtained with the 100 SA algorithm evolutions applying to simulated observations. Model denotes a likelihood model number, $\sigma_{\theta}$ denotes the empirical standard deviation, and $\Delta q$ denotes $q(0.975) - q(0.025)$.</caption><thead><tr><th>Parameter</th><th>Model</th><th>Mean ±2σ<sub>θ</sub></th><th>Δq</th></tr></thead><tbody><tr><td rowspan="4">a, km</td><td>1</td><td>10006.93 ± 197.86</td><td>710.27</td></tr><tr><td>2</td><td>9989.47 ± 97.75</td><td>341.65</td></tr><tr><td>3</td><td>9995.97 ± 132.85</td><td>516</td></tr><tr><td>4</td><td>9986.31 ± 93.55</td><td>380.68</td></tr><tr><td rowspan="4">e, –</td><td>1</td><td>0.5001 ± 0.02838</td><td>0.05194</td></tr><tr><td>2</td><td>0.4995 ± 0.00604</td><td>0.00992</td></tr><tr><td>3</td><td>0.4987 ± 0.02941</td><td>0.05947</td></tr><tr><td>4</td><td>0.4991 ± 0.00281</td><td>0.00562</td></tr><tr><td rowspan="4">i, <sup>a</sup> deg</td><td>1</td><td>134.98 ± 1.8</td><td>3.36</td></tr><tr><td>2</td><td>135.02 ± 0.08</td><td>0.15</td></tr><tr><td>3</td><td>135.04 ± 1.24</td><td>2.23</td></tr><tr><td>4</td><td>135.05 ± 0.33</td><td>0.65</td></tr><tr><td rowspan="4">Ω, <sup>a</sup> deg</td><td>1</td><td>44.95 ± 2.4</td><td>4.39</td></tr><tr><td>2</td><td>44.94 ± 1</td><td>2.1</td></tr><tr><td>3</td><td>44.94 ± 2.87</td><td>5.51</td></tr><tr><td>4</td><td>44.92 ± 1</td><td>2.54</td></tr><tr><td rowspan="4">ω, <sup>a</sup> deg</td><td>1</td><td>45 ± 2.01</td><td>3.64</td></tr><tr><td>2</td><td>45.04 ± 0.4</td><td>0.68</td></tr><tr><td>3</td><td>44.91 ± 1.3</td><td>2.74</td></tr><tr><td>4</td><td>45.07 ± 0.66</td><td>1.32</td></tr><tr><td rowspan="4">P, d</td><td>1</td><td>30 ± 0.52</td><td>0.93</td></tr><tr><td>2</td><td>29.99 ± 0.43</td><td>0.75</td></tr><tr><td>3</td><td>29.99 ± 0.44</td><td>0.93</td></tr><tr><td>4</td><td>29.99 ± 0.4</td><td>0.83</td></tr><tr><td rowspan="4">M, <sup>b</sup> deg</td><td>1</td><td>180 ± 7.93</td><td>14.42</td></tr><tr><td>2</td><td>179.86 ± 5.19</td><td>9.04</td></tr><tr><td>3</td><td>180.25 ± 5.12</td><td>9.6</td></tr><tr><td>4</td><td>179.78 ± 5.34</td><td>10.92</td></tr></tbody></table>

Notes. <sup>a</sup>Angles are referenced to the J2000 equatorial frame.

<sup>b</sup>The mean anomaly is referred to the epoch 15 d.

<table><caption>**Table 2.** Range of the average O – C residual, computed with equation (9), for obtained samples of orbits associated with the four likelihood models.</caption><thead><tr><th>Likelihood</th><th>min(O – C) – max(O – C), arcsec</th></tr></thead><tbody><tr><td>Model 1</td><td>0.00006–0.0156</td></tr><tr><td>Model 2</td><td>0.0001–0.0128</td></tr><tr><td>Model 3</td><td>0.0005–0.0126</td></tr><tr><td>Model 4</td><td>0.00005–0.0126</td></tr></tbody></table>

## 6.3 Real observations

For a real data case study, we chose the TNB Teharonhiawako (2001 QT<sub>297</sub>) with its companion called Sawiskera. The 16 observations published by Grundy et al. (2011) are used. The observations are distributed on the interval from 2001/10/11 to 2010/08/03, covering the entire orbital period (see Grundy et al. 2011). Some observations were made consecutively with a one day interval, but there are also those that follow after intervals of 1 yr and more (see Table A2). The estimated astrometric errors published with the observations were assumed to follow Gaussian distributions (Grundy et al. 2011) with standard deviation magnitudes from 0.13 to 0.003 arcsec.

First, the MAP algorithm for the aforementioned observations is used. Then, the LS orbit fitting is applied to the same data, using the initial set of parameters from Grundy et al. (2011). We consider the solution obtained with the LS method as a reference in order to verify the result from the MAP method. The reasons we use the LS techniques here are the following. The LS method is less computationally expensive. Although the LS method guarantees only local convergence, in the case of the good initial guess of parameters and the well-distributed data – here both given by Grundy et al. (2011) – the LS method provides an accurate solution, allowing us to consider it as a reference. Furthermore, the LS method allows us to easily switch between the not-weighted and weighted cases. Hence, an appropriate comparison can be made of the first/third and second/fourth MAP likelihood models with the not-weighted and weighted LS methods, respectively.

## Results

In the same way as for the simulated observations, for the real data, we obtain the distributions of each orbital parameter through 100 SA evaluations. For these distributions, we summarize the inter-quantile interval $q(0.025) – q(0.975)$ in Table 3. The set of parameters $\hat{\theta}$ of the best-fitting orbit with $2\sigma_{\theta}$ uncertainties and the LS result is also summarized in Table 3.

The results obtained for the real data using the SA algorithm are similar to those obtained for the simulated data. Again, the distribution shape, the numerical values of $\sigma_{\theta}$ and the length of the inter-quartile interval do not indicate an obvious bias, and the marginal distributions of the parameters tend to exhibit a Gaussian behaviour.

It can be noticed that the result, obtained with the first likelihood model, is slightly shifted from the LS result. The likelihood model 2 converges practically to the same result as the not-weighted LS. The slightly different results of models 1 and 2 (or LS) are not contradictory, but are explained by the different likelihood models. Concerning the WLS, the results have the same tendency as the not-weighted models: model 3 obtains resulting distributions shifted from the LS results and model 4 provides results very close to those of the LS (see Table 3).

Under the assumption that the LS method provides the true solution, the fact that the MAP algorithm gives results similar to the LS method validates the new method. Hereafter, we retain the second (for not-weighted observations) and the fourth (for weighted observations) likelihood models as more appropriate, because these models provide results closer to those of the LS method (in the not-weighted and weighted cases, respectively).

It should be noticed that for well-distributed observations and a good initial estimation of parameters, the LS method appears to be more practical, as it is less computationally expensive than the new proposed method. Hence, this particular situation does not require the use of our method; the classical LS-based method can be applied in order to save computation time. However, there is no rigorous criteria to state whether the data and initial guesses of parameters are good enough. Under this circumstance, the MAP method is more universal.

## 6.4 Ephemeris prediction

In order to verify the capability of the MAP method to make an ephemeris prediction, we apply it to a reduced set of observations. Namely, the eight first observations from the 16 given in Grundy
<table><caption>**Table 3.** Derived parameters $\hat{\theta}$ with uncertainty estimate for MAP and LS methods obtained from real observations of Teharonhiawako 2001 QT$_{297}$. Model denotes a likelihood model of the MAP method or LS and WLS methods, $\sigma_{\theta}$ denotes an empirical standard deviation and $\Delta q$ denotes $q(0.975) - q(0.025)$.</caption><thead><tr><th>Parameter</th><th>Model</th><th>θ̂ ± 2 σ<sub>θ</sub></th><th>Δq</th></tr></thead><tbody><tr><td rowspan="5">a, km</td><td>1</td><td>28100.87 ± 90.06</td><td>180.11</td></tr><tr><td>2</td><td>28124.49 ± 1.69</td><td>3.37</td></tr><tr><td>3</td><td>28006.11 ± 143.91</td><td>287.82</td></tr><tr><td>4</td><td>27777.15 ± 2.4</td><td>4.79</td></tr><tr><td>LS</td><td>28125.8 ± 407.2</td><td></td></tr><tr><td></td><td>WLS</td><td>27780.18 ± 1.78</td><td></td></tr><tr><td rowspan="5">e, -</td><td>1</td><td>0.2439 ± 0.0017</td><td>0.003</td></tr><tr><td>2</td><td>0.2435 ± 0.00006</td><td>0.0001</td></tr><tr><td>3</td><td>0.249 ± 0.00213</td><td>0.00373</td></tr><tr><td>4</td><td>0.2547 ± 0.00014</td><td>0.00027</td></tr><tr><td>LS</td><td>0.2436 ± 0.0081</td><td></td></tr><tr><td></td><td>WLS</td><td>0.2548 ± 4e-5</td><td></td></tr><tr><td rowspan="5">i, <sup>a</sup> deg</td><td>1</td><td>145.49 ± 0.26</td><td>0.45</td></tr><tr><td>2</td><td>144.01 ± 0.01</td><td>0.02</td></tr><tr><td>3</td><td>145.37 ± 0.44</td><td>0.74</td></tr><tr><td>4</td><td>143.98 ± 0</td><td>0.01</td></tr><tr><td>LS</td><td>144.01 ± 1.5</td><td></td></tr><tr><td></td><td>WLS</td><td>143.99 ± 0.01</td><td></td></tr><tr><td rowspan="5">Ω, <sup>a</sup> deg</td><td>1</td><td>54.28 ± 0.71</td><td>1.3</td></tr><tr><td>2</td><td>51.88 ± 0.02</td><td>0.04</td></tr><tr><td>3</td><td>55.11 ± 0.42</td><td>0.77</td></tr><tr><td>4</td><td>55.08 ± 0.02</td><td>0.07</td></tr><tr><td>LS</td><td>51.89 ± 2.5</td><td></td></tr><tr><td></td><td>WLS</td><td>55.09 ± 0.01</td><td></td></tr><tr><td rowspan="5">ω, <sup>a</sup> deg</td><td>1</td><td>324.97 ± 0.71</td><td>1.3</td></tr><tr><td>2</td><td>324.32 ± 0.03</td><td>0.05</td></tr><tr><td>3</td><td>325.65 ± 0.2</td><td>0.39</td></tr><tr><td>4</td><td>324.81 ± 0.02</td><td>0.04</td></tr><tr><td>LS</td><td>324.34 ± 2.5</td><td></td></tr><tr><td></td><td>WLS</td><td>324.84 ± 0.01</td><td></td></tr><tr><td rowspan="5">P, d</td><td>1</td><td>828.21 ± 0.36</td><td>0.7</td></tr><tr><td>2</td><td>828.17 ± 0.01</td><td>0.03</td></tr><tr><td>3</td><td>828.4 ± 0.24</td><td>0.39</td></tr><tr><td>4</td><td>828.07 ± 0.02</td><td>0.04</td></tr><tr><td>LS</td><td>828.15 ± 80.9</td><td></td></tr><tr><td></td><td>WLS</td><td>828.07 ± 0.004</td><td></td></tr><tr><td rowspan="5">M, <sup>b</sup> deg</td><td>1</td><td>276.02 ± 0.24</td><td>0.46</td></tr><tr><td>2</td><td>275.27 ± 0.02</td><td>0.03</td></tr><tr><td>3</td><td>276 ± 0.18</td><td>0.38</td></tr><tr><td>4</td><td>276.7 ± 0.02</td><td>0.03</td></tr><tr><td>LS</td><td>275.26 ± 1.5</td><td></td></tr><tr><td></td><td>WLS</td><td>276.68 ± 0.01</td><td></td></tr></tbody></table>

Notes. <sup>a</sup>Angles are referenced to the J2000 equatorial coordinates.

<sup>b</sup>The mean anomaly is referred to the epoch JD 245 2000.0.

et al. (2011) are used for the orbit determination problem. It should be noted that for this reduced set of observations, the LS method does not converge, whereas the MAP method provides a solution. Using this solution, we calculate the ephemeris at the epochs of the other eight observations following after the first set. Following the method presented in Section 5.5, 95 per cent inter-quantile intervals were computed for the calculated positions.

Figs 4 and 5, associated with the second and fourth likelihood models, respectively, show observed and predicted positions to-

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Ephemeris prediction for the result of the MAP method with the likelihood model 2. The calculated positions (black circles) are compared with given observed positions (crosses) by the x and y coordinates on the sky plane. Black bars denote the interval of 2.5–97.5 per cent quantiles. The dotted line corresponds to the calculated positions for the orbit, obtained with the entire set of observations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Ephemeris prediction for the result of the MAP method with the likelihood model 4. The calculated positions (black circles) are compared with given observed positions (crosses) by the x and y coordinates on the sky plane. Black bars denote the interval of 2.5–97.5 per cent quantiles. The dotted line corresponds to the calculated positions for the orbit, obtained with the entire set of observations.</figcaption></figure>

gether with the 95 per cent inter-quantile interval. It can be seen, for both likelihood models, that the estimated uncertainties can be as large as the apparent size of orbit (see the x and y amplitudes in the figures). The test of the validity of the model is rejected only by one observation: the ninth (see the y coordinate of the first observed
position in Figs 4 and 5), which can be caused by an observational bias or the imperfection of the chosen physical model. Still, all of the other observed positions are situated within the range of the estimated interval. Hence, we can consider the observed positions as realizations of the stochastic process providing the different MAP solutions.

## 7 CONCLUSION AND PERSPECTIVES

The estimation of orbital parameters using the MAP method with the SA algorithm has been proposed for the binary asteroid orbit determination problem. The new method has been successfully applied to simulated and real observations and has been verified for ephemeris prediction. The developed algorithm is prepared to be implemented into the Gaia mission (Prusti et al. 2016) data-reduction pipeline.

The MAP method is independent of the initial orbit estimation, guaranteeing the theoretical convergence towards the global optimum solution. This advantage allows us to use it for newly discovered asteroids, when an initial parameter estimation is difficult. At the same time, if any available information helps to constrain the parameter search space, it can be easily introduced through the a priori term.

The new method enables us to consider and to compare different forms of observational error distribution, implementing it by the likelihood modelling. This aspect can be particularly useful when the observational errors cannot be assumed to obey Gaussian distribution.

The sampling mechanism within the SA algorithm is done through the orbital parameters. This approach simplifies the entire procedure, as no methods (such as the Thiele-Innes method, or equivalent) are required for computation of orbital parameters from observations.

The parameter estimation by the MAP method with the SA algorithm implementation is still developing. The method has been implemented for only a few binaries (see also Kovalenko et al. (2016)) and we expect to improve it when applied in more various cases. For this purpose, the following perspectives are depicted. First, even if guaranteeing the convergence of the method and the quality of the obtained results, the choices of the cooling schedule and of the proposal distribution in the SA algorithm are not optimal. Second, the dynamical model of a binary motion can be completed by considering a perturbed orbit instead of Keplerian. Third, the new method can be expanded for the non-resolved binary asteroids.

## ACKNOWLEDGEMENTS

This work is supported by Labex ESEP (ANR N 2011-LABX-030) and by the Russian Foundation for Basic Research, project no. 16-52-150005-CNRS-a. The work of RSS was carried out while he was working at Université de Lille – Laboratoire Paul Painlevé.

## REFERENCES

* Aitken R. G., 1964, The Binary Stars. Dover, New York

* Descamps P., 2005, Celestial Mechanics and Dynamical Astronomy, 92, 381

* Emelyanov N., 2005, A&A, 435, 1173

* Geman S., Geman D., 1984, IEEE Transactions on Pattern Analysis and Machine Intelligence, 6, 721

* Grundy W. et al., 2008, Icarus, 197, 260

* Grundy W. M. et al., 2011, Icarus, 213, 678

* Hastings W. K., 1970, Biometrika, 57, 97

* Hestroffer D., Vachier F., Balat B., 2005, Earth, Moon, and Planets, 97, 245

* Jeffreys H., 1946, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 186, 453

* Jeffreys H., 1973, Scientific Inference. Cambridge University Press, Cambridge

* Johnston W. R., 2014, NASA Planetary Data System, 219

* Kirkpatrick S., Gelatt C. D., Vecchi M., 1983, Science, 220, 671

* Kovalenko I. D., Stoica R. S., Emelyanov N. V., Doressoundiram A., Hestroffer D., 2016, in Chesley S. R., Morbidelli A., Jedicke R., Farnocchia D., eds, Proc. IAU Symp. Vol. 318, Asteroids: New Observations, New Models. Kluwer, Dordrecht, p. 244

* Li S. Z., 2009, Markov Random Field Modelling in Image Analysis. Springer Science & Business Media, Berlin

* Muinonen K., Powell E., 1993, Icarus, 104, 255

* Muinonen K., Granvik M., Oszkiewicz D., Pieniluoma T., Pentikäinen H., 2012, Planetary and Space Science, 73, 15

* Nelder J. A., Mead R., 1965, The Computer Journal, 7, 308

* Noll K. S., Stephens D. C., Grundy W. M., Griffin I., 2004, Icarus, 172, 402

* Oszkiewicz D., Muinonen K., Virtanen J., Granvik M., 2009, Meteoritics and Planetary Science, 44, 1897

* Oszkiewicz D., Hestroffer D., David P., 2013, in European Planetary Science Congress 2013 (information online at http://meetings.copernicus.org/epsc2013), id. EPSC2013-522

* Prusti T. et al., 2016, A&A, 595, A1

* Roberts G. O., Smith A. F. M., 1994, Stochastic Processes and their Applications, 49, 207

* Roberts G. O., Tweedie R. L., 1996, Biometrika, 83, 95

* Stoica R. S., Gregori P., Mateu J., 2005, Stochastic Processes and their Applications, 115, 1860

* Vachier F., Berthier J., Marchis F., 2012, A&A, 543, A68

* Virtanen J., Muinonen K., Powell E., 2001, Icarus, 154, 412

## APPENDIX: SIMULATED AND REAL OBSERVATIONS

<table><caption>**Table A1.** Simulated observations: $x = (\alpha_2 - \alpha_1)\cos \delta_1$ and $y = \delta_2 - \delta_1$, where $\alpha$ is right ascension and $\delta$ is declination of the target asteroid, referring to primary (1) and secondary (2), respectively, and $R$ is the distance from the observer to the target. Estimated $\sigma_{x,i}$ and $\sigma_{y,i}$ uncertainties in the final digits are indicated in parentheses.</caption><thead><tr><th>Day</th><th>x (arcsec)</th><th>y (arcsec)</th><th>R (au)</th><th>α (deg)</th><th>δ (deg)</th></tr></thead><tbody><tr><td>5</td><td>-0.0072 (10)</td><td>0.1741 (10)</td><td>44.87</td><td>56.02</td><td>24.01</td></tr><tr><td>14</td><td>0.2663 (20)</td><td>-0.0662 (10)</td><td>44.74</td><td>56.18</td><td>24.01</td></tr><tr><td>26</td><td>0.1413 (30)</td><td>-0.2143 (40)</td><td>44.60</td><td>56.02</td><td>24.01</td></tr><tr><td>30</td><td>-0.0939 (40)</td><td>0.0318 (10)</td><td>44.56</td><td>56.21</td><td>24.01</td></tr><tr><td>38</td><td>0.1068 (50)</td><td>0.1118 (40)</td><td>44.49</td><td>56.07</td><td>24.01</td></tr><tr><td>41</td><td>0.2004 (90)</td><td>0.0248 (10)</td><td>44.46</td><td>56.02</td><td>24.00</td></tr><tr><td>47</td><td>0.3012 (10)</td><td>-0.1494 (40)</td><td>44.43</td><td>55.15</td><td>24.02</td></tr><tr><td>50</td><td>0.2989 (70)</td><td>-0.2142 (10)</td><td>44.41</td><td>55.09</td><td>24.01</td></tr><tr><td>59</td><td>-0.0405 (10)</td><td>-0.0545 (40)</td><td>44.39</td><td>55.16</td><td>24.01</td></tr><tr><td>62</td><td>-0.1082 (80)</td><td>0.1558 (10)</td><td>44.38</td><td>55.10</td><td>24.01</td></tr></tbody></table>
<table><caption>**Table A2.** Observations of (2001 QT297) Teharonhiawako (Grundy et al. 2011), with the relative $x = (\alpha_2 - \alpha_1)\cos\delta_1$ and $y = \delta_2 - \delta_1$, where $\alpha$ and $\delta$ are right ascension and declination respectively, and subscripts 1 and 2 refer to the primary and secondary, respectively. The distance from the observer to the target is $R$. Estimated $\sigma_{x,i}$ and $\sigma_{y,i}$ uncertainties in the final digits are indicated in parentheses.</caption><thead><tr><th>UT day</th><th>UT hour</th><th>R (au)</th><th>x (arcsec)</th><th>y (arcsec)</th></tr></thead><tbody><tr><td>2001/10/11</td><td>0.9528</td><td>44.370</td><td>+0.5390 (51)</td><td>-0.2770 (52)</td></tr><tr><td>2001/10/12</td><td>1.8730</td><td>44.385</td><td>+0.5460 (81)</td><td>-0.2675 (84)</td></tr><tr><td>2001/11/01</td><td>0.7697</td><td>44.699</td><td>+0.624 (21)</td><td>-0.214 (24)</td></tr><tr><td>2001/11/02</td><td>0.4299</td><td>44.716</td><td>+0.644 (21)</td><td>-0.184 (26)</td></tr><tr><td>2001/11/03</td><td>0.3249</td><td>44.733</td><td>+0.642 (40)</td><td>-0.193 (39)</td></tr><tr><td>2001/11/04</td><td>0.9046</td><td>44.750</td><td>+0.645 (21)</td><td>-0.138 (35)</td></tr><tr><td>2002/07/13</td><td>6.7387</td><td>44.130</td><td>-0.314 (23)</td><td>+0.692 (29)</td></tr><tr><td>2002/07/18</td><td>6.9538</td><td>44.085</td><td>-0.344 (68)</td><td>+0.700 (55)</td></tr><tr><td>2002/08/07</td><td>4.5629</td><td>43.970</td><td>-0.43 (13)</td><td>+0.81 (13)</td></tr><tr><td>2002/09/08</td><td>5.7632</td><td>44.024</td><td>-0.658 (91)</td><td>+0.658 (91)</td></tr><tr><td>2003/10/23</td><td>1.7567</td><td>44.536</td><td>-0.012 (60)</td><td>-0.527 (50)</td></tr><tr><td>2004/05/25</td><td>8.7890</td><td>44.878</td><td>+0.4350 (70)</td><td>+0.4560 (70)</td></tr><tr><td>2004/09/13</td><td>3.2531</td><td>44.074</td><td>-0.1330 (70)</td><td>+0.6990 (60)</td></tr><tr><td>2005/07/11</td><td>5.8782</td><td>44.234</td><td>-1.0020 (80)</td><td>-0.0440 (80)</td></tr><tr><td>2009/12/12</td><td>5.2070</td><td>45.368</td><td>-1.0257 (77)</td><td>+0.1098 (53)</td></tr><tr><td>2010/08/03</td><td>10.1942</td><td>44.166</td><td>-0.0032 (30)</td><td>-0.5015 (30)</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.