Romain Serra
email: serra@laas.fr

Denis Arzelier
email: arzelier@laas.fr

Mioara Joldes
email: mjoldes@laas.fr

Jean-Bernard Lasserre
email: lasserre@laas.fr

Aude Rondepierre
email: aude.rondepierre@math.univ-toulouse.fr

Bruno Salvy
email: bruno.salvy@inria.fr

A Fast and Accurate Method to Compute the Probability of Collision for Short-term Space Encounters

This article provides a new method for computing the probability of collision between two spherical space objects involved in a short-term encounter under Gaussiandistributed uncertainty. In this model of conjunction, classical assumptions reduce the probability of collision to the integral of a two-dimensional Gaussian probability density function over a disk. The computational method presented here is based on an analytic expression for the integral, derived by use of Laplace transform and D-finite functions properties. The formula has the form of a product between an exponential term and a convergent power series with positive coefficients. Analytic bounds on the truncation error are also derived and are used to obtain a very accurate algorithm.

orbiting spacecraft. Avoiding collision has become a usual and necessary task for many satellites on duty. In order to prevent collisions, space debris are nowadays radar-tracked and conjunctions with operational spacecraft can be predicted. From this information on a possible collision between two objects, it is then the task of the operator to assess the risk. Because of the uncertain nature of the data, the decision parameter is often computed as a probability of collision. If it exceeds some tolerance threshold, an evasive maneuver will be performed. When modeling conjunctions, two assumptions are usually made in the literature [START_REF] Alfriend | Probability of collision error analysis[END_REF][START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF]. First, the geometrical shape of the two objects is supposed to be spherical. Second, uncertainty on the state vectors of the two objects is modeled with independent normal probability distributions.

Moreover, in this article, only short-term encounters are considered i.e., conjunctions with high relative velocities. Such cases typically arise in Low Earth Orbit where orbital velocities are fairly high. Two assumptions define the short-term encounter model [START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF]. First, the relative trajectory of the objects during the encounter is considered to be a straight line [START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF]. Secondly, uncertainty on the velocities is neglected. The probability of collision can then be expressed as a two-dimensional Gaussian integral over a disk. The integrand is the probability density function of the relative coordinates in the so-called encounter plane [START_REF] Chan | Spacecraft Collision Probability[END_REF] while the domain of integration is the cross-section of the spherical combined object. As a result, the only inputs needed for the computation of the probability of collision are: the radius of the combined spherical object, the miss distance and the covariance matrix of the relative position in the encounter plane at reference time. A more detailed description is given in Section II.

The computation of the short-term encounter probability of collision has been dealt with in several ways. First, as a probability, it can be obtained from Monte Carlo trials [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF], but this is a time-consuming process. Another class of methods are based on numerical integration schemes.

In chronological order of apparition, they are: Foster's [START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF], Patera's [START_REF] Patera | General method for calculating satellite conjunction probability[END_REF] and Alfano's [START_REF] Alfano | Aerospace support to space situation awareness[END_REF]. The main disadvantage of this family of methods is that they are strongly dependent upon the chosen integration method and suppose to manage a sensitive trade-off between precision and computation time. Finally, the approach closest to ours is Chan's [START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Chan | Collision probability analyses for earth-orbiting satellites[END_REF] who derives an analytic formula, but with an approximation with respect to the initial model.

Contribution

In Section III a new analytical formula is proposed, with the following distinguishing features:

• It has the form of a product between an exponential term and a convergent power series with positive terms. The series is obtained by use of Laplace transform techniques [START_REF] Lasserre | Solving a class of multivariate integration problems via laplace techniques[END_REF]. The explicit form of the linear recurrence satisfied by the coefficients of the series is derived using properties of D-finite functions i.e., solutions of linear ordinary differential equations with polynomial coefficients [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF][START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF].

• Analytical bounds on the truncation error are provided, enabling the user to evaluate the truncated series as accurately as required. As a by-product of these bounds on the truncation error, one obtains analytical bounds on the probability of collision.

• In contrast to Chan's method, it does not make any approximation of the integral to be computed. In is shown in Section IV that Chan's formula is a special case of the new formula.

• When compared to other methods on a large number of test cases, the resulting algorithm performs well in practice. As illustrated in Section V by numerical examples, it allows for a very accurate and fast way to estimate the risk in most cases, which is highly desirable from an operator point of view.

II. Encounter Model

In this section, the general framework of on-orbit collision risk is briefly described. Next, we present the model of short-term encounters between two spherical objects under Gaussiandistributed uncertainty.

A. General encounter between two space objects

The system considered is composed of two space objects subject to orbital dynamics over some time interval. One is called the primary object p and the other one the secondary object s. Typically, the primary is an active spacecraft while the secondary is a space debris. Information -such as radar measurements -describes their respective geometry, position and velocity at some detection time.

The initial conditions -position and velocity of p and s at reference time -are supposed to be uncertain parameters, resulting in uncertain dynamics for each object. The reference time usually belongs to the time interval when both objects are in the collision range. Generally, it is the so-called nominal time of closest approach (TCA) [START_REF] Akella | Probability of collision between space objects[END_REF]. This date is defined as the instant when the mean relative distance between the two objects reaches its global minimum over time. This minimum value is another quantity of great use referred to as the nominal miss distance [START_REF] Akella | Probability of collision between space objects[END_REF]. Considering that statistical assumptions are usually made on these initial conditions, the collision risk is naturally characterized via a probabilistic metric. As a result, it can be defined as the integral of a probability density function (pdf) over some domain of integration. Without any further assumption, computing the probability of collision is a hard integration problem for two main reasons. First, the domain of integration is often impossible to derive analytically, because of its dependency upon the dynamics of each object which may be fairly complex depending on the model being considered. Secondly, integrating the probability density function over this domain is also difficult in general since, even for instance with Gaussian density functions, such an integral generally admits no closed-form formula [START_REF] Genz | Computation of multivariate normal and t probabilities[END_REF].

A general method to handle the computation of the probability of collision is based on Monte-Carlo trials [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF]. It consists in simulating the dynamics of a large number of scenarios and in counting the number of effective collisions. Unfortunately, the number of required trials goes up with the targeted precision and makes this method quite time-consuming. For so-called rare events -events with very low probability values -it can even become intractable [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF]. Therefore, simplified encounter models have been designed to reduce the complexity of the computation of the probability of collision, allowing for faster dedicated algorithms.

B. Short-term encounter between spherical objects

One of these simplified encounter models is the so-called short-term encounter model [START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF][START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF] that can be used to describe a number of conjunctions when the magnitude of relative velocities at stake is sufficiently high. This typically happens in Low-Earth Orbit (LEO) where orbital velocities are the largest in norm [START_REF] Chan | Spacecraft Collision Probability[END_REF]. Investigations about the range of validity of the short-term encounter model can be found in [START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Evaluating the Short Encounter Assumption of the Probability of Collision Formula[END_REF]. In the case of a short-term encounter between spherical objects under Gaussian-distributed uncertainty, the probability of collision reduces to a 2-D Gaussian integral over a disk [START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF][START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF]. The present article addresses the practical computation of this quantity.

Five assumptions are needed to define the short-term encounter model under Gaussian-distributed uncertainty for spherical objects:

1. The relative trajectories are approximated as rectilinear.

2. The velocities at reference time are considered as deterministic variables.

3. Initial position vectors of both objects are Gaussian independent random vectors.

4. Each object is approximated by a spherical geometrical shape.

5. The time boundaries of the conjunction are extended to infinity.

Encounter frame

The fact that the relative motion is rectilinear motivates the choice of a frame of study with one axis along the relative velocity. One possibility is to introduce the so-called encounter frame [START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF] defined at reference time. Its exact definition depends on the authors (see for instance [START_REF] Akella | Probability of collision between space objects[END_REF] and [START_REF] Chan | Spacecraft Collision Probability[END_REF] for two different frames) but it always has the two following characteristics: it is centered on the mean position of one of the two objects and is built from the so-called encounter plane. This plane contains the origin and is orthogonal to the direction of the relative velocity.

The configuration considered in this study is represented in Figure 1. The origin of the frame is located at the center of the primary object p. The basis vector e z is oriented along the relative velocity v = v s -v p . The basis vector e x belongs to the encounter plane: it points towards the orthogonal projection of the mean relative position µ(r 0 ) onto the encounter plane. Finally, the basis vector e ỹ completes the right-handed system and thus belongs to the encounter plane as well.

In summary, one has

e z = v v , e ỹ = v × µ(r 0 ) v × µ(r 0 ) , e x = e ỹ × e z . (1)
Let (x m , 0, zm ) be the coordinates of the mean relative position in the encounter frame where the ỹ coordinate is zero by construction. It is worth noticing that, since the relative trajectory is rectilinear, xm is in fact equal to the miss distance.

Integral representation

Under the assumptions of the short-term encounter model, the probability of collision can be formulated as a 2-D integral in the encounter plane. For spherical objects, the domain of integration is a closed disk B((0, 0), R) centered at the origin of radius R. The quantity R is the combined radius and is defined as the sum of the respective radii of the two objects i.e. R = R p + R s .

The two-dimensional pdf involved in the probability of collision describes the distribution of the relative position in the encounter plane. From the hypothesis on the nature of uncertainty, it is a multivariate normal law. Therefore, it is completely defined by its mean vector and its variancecovariance matrix [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF]. Let Σ xỹ be the covariance matrix of the relative coordinates in the encounter plane. The probability of collision can then be written as:

P c = 1 2π |Σ xỹ | B((0,0),R) exp - 1 2 [x -xm ỹ] Σ -1 xỹ [x -xm ỹ] T dxdỹ. (2)
Equation [START_REF] Akella | Probability of collision between space objects[END_REF] shows that the probability of collision only depends on the combined radius R, the miss distance xm and the covariance matrix Σ xỹ of the relative coordinates in the encounter plane.

Frame rotation

Further calculations allow to write a formula with a simpler integrand. As a covariance matrix, Σ xỹ can be written as:

Σ xỹ =     σ 2 x ρ xỹ σ xσ ỹ ρ xỹ σ xσ ỹ σ 2 ỹ     , (3)
where σ x, σ ỹ ∈ R + * are the standard deviations of the relative coordinates in the encounter plane and

ρ xỹ ∈ (-1, 1
) is the correlation coefficient. In order to eliminate the cross-terms of the Gaussian function, a rotation to the principal axis of the covariance matrix is performed in the encounter plane (see Figure 2). The new coordinates, denoted (x, y), are respectively along the major and the minor axis. Consistently with Chan's notation [START_REF] Chan | Spacecraft Collision Probability[END_REF], let -θ be the corresponding rotation angle.

Nevertheless, the reader should be aware that Chan actually uses this symbol to denote any angle rotating to the principal axis, not necessarily with the major axis aligned with the first axis as done here. This transformation does not change the nature of the domain of integration which remains a disk of radius R centered at the origin. It follows that:

P c = 1 2πσ x σ y B((0,0),R) exp - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dxdy, (4)
where the quantities σ x and σ y are standard deviations of the new coordinates. As a matter of fact, σ 2

x and σ 2 y are respectively the largest and the smallest eigenvalues of Σ xỹ :

     σ 2 x + σ 2 ỹ 2 ± σ 2 x -σ 2 ỹ 2 2 + ρ 2 xỹ σ 2 xσ 2 ỹ      . ( 5
)
Accordingly to the rotation of angle -θ, one has:

x m = xm cos θ, y m = -x m sin θ. ( 6
)
The angle θ can always be chosen in (-π 2 , π 2 ]. If ρ xỹ = 0, it is given by the following formula:

θ = arctan    σ 2 ỹ -σ 2 x 2ρ xỹ σ xσ ỹ + sign(ρ xỹ ) 1 + σ 2 ỹ -σ 2 x 2ρ xỹ σ xσ ỹ 2    . ( 7
)
If ρ xỹ = 0, the matrix Σ xỹ is in fact already diagonal and two subcases can occur. If σ x ≥ σ ỹ , no rotation is needed and θ = 0; if σ x < σ ỹ , one only needs to interchange the axis and θ = π 2 .

x The contribution of the article is to propose a new systematic computation method that will be presented in details in the next section.

C. State of the art for the computation of Pc

Monte Carlo methods can be used to compute the probability of collision for short-term encounters. But as for the general collision model, they are time-consuming and ill-fitted to compute probabilities with very low values. Moreover, they do not take advantage of the particular expression of the short-term encounter formula. On the other hand, four methods have been developed so far to specifically handle this computation.

The first of these methods -used by NASA [START_REF] Chan | Spacecraft Collision Probability[END_REF] -originates from an article by Foster and Estes [START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF].

It consists in a direct application of a numerical discretization scheme to the two-dimensional integral after switching to polar coordinates. While being fairly precise, it is relatively slow compared to other methods [START_REF] Chan | Spacecraft Collision Probability[END_REF]. Patera [START_REF] Patera | General method for calculating satellite conjunction probability[END_REF] reformulates the probability of collision as a one-dimensional integral by writing it as a path-integral over the contour of the domain of integration. A numerical scheme -namely the midpoint method -is used for actual computation. By use of the error function,

Alfano [START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF] also formulates the probability of collision as a one-dimensional integral. Then, Simpson's integration scheme is applied to compute it. The outcome is a finite sum involving exponential and error functions.

Chan's method [START_REF] Chan | Spacecraft Collision Probability[END_REF] is the only one based on an analytical formula rather than a numerical scheme. However, it is derived via an approximation of the domain of integration. Geometrically, Chan considers the integration problem of an isotropic Gaussian function over an ellipse -isotropic meaning that the level sets of the Gaussian function are circles. It is strictly equivalent to the original integration problem of an anisotropic Gaussian function over a disk. Then, the actual elliptic domain of integration is approximated by a disk of equivalent area and same center. Thus, the probability of collision reduces to the integral of an isotropic Gaussian function over a disk shifted from its peak. This integration problem can be solved analytically by use of properties of the first modified Bessel function. Finally, Chan is able to write the probability of collision as an infinite series with positive terms. The advantage of this approach is to provide the user with a truncation error that can be bounded. The major drawback lies in the fact that the formula itself is only an approximation of the original integral and that the discrepancy between the two is hard to quantify. One of the advantages of the present paper is to propose a different analytic formula for the probability of collision -derived without any approximation.

III. A new method to compute the probability of collision

The three main ingredients used in the derivation of the new formula are:

1. A method introduced by Lasserre and Zeron [START_REF] Lasserre | Solving a class of multivariate integration problems via laplace techniques[END_REF] to integrate Gaussian functions over Euclidean balls is used. Specifically, equation ( 4) is reformulated as a function g(ξ), with ξ = R 2 , whose Laplace transform L g is computed in closed-form. The function is then expanded in a Laurent series. A term by term application of the inverse Laplace transform leads to a power series for the initial integral.

2. For the first time, a simple form for these coefficients is found. In Section III B, they are proved to satisfy a linear recurrence with polynomial coefficients (in the index variable). Classical theory and properties of D-finite (or holonomic) functions [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF][START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF] are used for this purpose.

3. From a numerical point of view, the direct evaluation of the series obtained for g(ξ) can be difficult. Roughly speaking, although the power series expansion of g(ξ) is convergent, the evaluation of the sum in finite precision arithmetic is prone to high cancellation [START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[END_REF][START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF]. This comes from consecutive terms that are close in magnitude, but of different signs, so that their sum in finite precision arithmetic contains very few correct significant digits. This makes the power series evaluation impractical for large values of ξ. Therefore, a so-called preconditionning [START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[END_REF][START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF] is used in order to obtain a series suitable for low cancellation numerical evaluation. Instead of g, the function ψ • g is considered, where the so-called preconditionner ψ is a D-finite function as well. By carefully choosing the preconditionner, both ψ and ψ • g can be efficiently numerically evaluated.

In the end, the formula for the probability of collision has the following form:

P c = exp - R 2 2σ 2 y +∞ k=0 β k R 2k , (8)
where (β k ) k≥0 is a positive sequence given by an explicit linear recurrence. This is detailed in what follows.

A. Laplace transform computation of the preconditioned probability of collision

The probability of collision given by (4) may be written as P c = g(R 2 ) where the function g : R + → R + is defined as follows:

g(ξ) = 1 2πσ x σ y B((0,0), √ ξ) exp - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dxdy. ( 9
)
This function is in fact the cumulative density function (cdf) of the random variable Ξ = X 2 + Y 2

where X ∼ N (x m , σ 2 x ) and Y ∼ N (y m , σ 2 y ) are independent random variables.

The function f = ψ • g is considered instead of the function g to prevent the cancellation phenomenon. The heuristic for a good choice of ψ is based on a method presented in [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF] (using several complex analysis results and properties obtained for g). In our case, it may be proved that an exponential function ψ : ξ → exp(pξ) is a good choice for several values of p. The perconditioner

p = 1 2σ 2 y
is chosen, since it leads to a simple proof of positivity of the series coefficients of g and moreover, allows for a unified and generalized way of interpreting the exact analytic formula obtained by Chan [START_REF] Chan | Spacecraft Collision Probability[END_REF]Chap. 4] in the isotropic case i.e., for σ x = σ y (see Section IV). The following proposition provides the closed-form Laplace transform for f .

Proposition III.1. Let f : R + → R + , such that

f (ξ) = exp(pξ)g(ξ), ( 10
)
where p ∈ R + . Then, the Laplace transform L f of the function f , is given by

L f (λ) = exp -(λ -p) x 2 m 2(λ-p)σ 2 x +1 + y 2 m 2(λ-p)σ 2 y +1 (λ -p) 2(λ -p)σ 2 x + 1 2(λ -p)σ 2 y + 1 , for |λ| > p. ( 11
)
Proof. The function g can be expressed as an integral over R 2 by use of the set-indicator function:

g(ξ) = 1 2πσ x σ y R 2 1 B((0,0), √ ξ) (x, y) exp - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dxdy. (12)
Consider now the Laplace transform L g of g. Then, for all λ ∈ C such that Re(λ) > 0,

L g (λ) = +∞ 0 g(ξ) exp(-λξ)dξ, = 1 2πσ x σ y +∞ 0 R 2 1 B((0,0), √ ξ) (x, y) exp -λξ - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dxdydξ. ( 13
)
Interchanging integrals is justified by Fubini's theorem, by absolute convergence of [START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF], that becomes:

L g (λ) = 1 2πσ x σ y R 2 +∞ 0 1 B((0,0), √ ξ) (x, y) exp (-λξ) exp - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dξ dxdy, = 1 2πσ x σ y λ R 2 exp -λ(x 2 + y 2 ) - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y dxdy, ( 14
) = 1 2πσ x σ y λ +∞ -∞ exp -λx 2 - (x -x m ) 2 2σ 2 x dx +∞ -∞ exp -λy 2 - (y -y m ) 2 2σ 2 y dy. ( 15
)
Next, each individual one-dimensional integral can be computed analytically. Indeed, after a change of variables aiming at completing the square and scaling the expression inside the exponential, the

remaining integral is +∞ -∞ exp - r 2 2 dr = √ 2π.
In the end, one gets:

L g (λ) = exp -λ x 2 m 2λσ 2 x +1 + y 2 m 2λσ 2 y +1 λ 2λσ 2 x + 1 2λσ 2 y + 1 . ( 16
)
Multiplying the function g(ξ) by exp(pξ) corresponds to a translation in the Laplace domain

L f (λ) = L g (λ -p).
Plugging this p-shift into ( 16) gives [START_REF] Lasserre | Solving a class of multivariate integration problems via laplace techniques[END_REF].

Note that in fact the function λ → λL g (-λ) is, by definition, the moment generating function of the random variable Ξ previously defined [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF].

B. Power series of f using D-finiteness

For simplicity, the following notations are used:

p = 1 2σ 2 y , φ = 1 - σ 2 y σ 2 x , ω x = x 2 m 4σ 4 x , ω y = y 2 m 4σ 4
y and α 0 = 1 2σ x σ y exp - 1 2 x 2 m σ 2 x + y 2 m σ 2 y .
Note that 0 ≤ φ < 1, ω x ≥ 0, ω y ≥ 0 and α 0 > 0. With this notation, equation [START_REF] Lasserre | Solving a class of multivariate integration problems via laplace techniques[END_REF] becomes:

L f (λ) = α 0 exp ωy λ + ωx λ-pφ λ(λ -pφ)(λ -p) . ( 17
)
From [START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[END_REF], it is possible to obtain the power series of f by performing a term-by-term inverse Laplace Transform on the principal part of Laurent series (i.e., the series of terms with negative degree) of L f [19, Chap. 9], [20, Chap. 2.14].

Note first that the principal part of the Laurent series of L f can be obtained by computing the Taylor expansion at infinity of L f . This is done by computing the coefficients of the power series expansion (at zero) of L f (λ -1 ). Since the first two coefficients of L f (λ -1 ) are zero, let Lf (λ) := λ -2 L f (λ -1 ) in order to simplify subsequent formulas. The multiplication by λ -2 is a technicality: it means just a shift by two positions in the series coefficients.

A simple computation gives

L f (λ -1 ) = α0λ 2 exp ωyλ-ωx pφ - ωx pφ(pφλ-1) √ 1-pφλ(1-pλ)
. A simple derivation shows that Lf verifies:

d Lf (λ) dλ = ϕ(λ) Lf (λ), Lf (0) = α 0 , ( 18
)
where

ϕ(λ) = ω y - pφ 2(pφλ -1) - p pλ -1 + ω x (pφλ -1) 2 . ( 19
)
Hence, Lf (λ) is D-finite and may be expanded in ∞ k=0 α k λ k , for all {λ, |λ| < 1 p }. Therefore, the sequence (α k ) k≥0 satisfies a linear recurrence obtained using Gfun and given by:

(k + 4)α k+4 = -p 3 φ 2 ω y α k + p 2 φ pφ k + 5 2 + 2ω y φ 2 + 1 α k+1 -p pφ φ 2 + 1 (2k + 5) + φ 2ω y + 3p 2 + ω x + ω y α k+2 + p(2φ + 1)(k + 3) + p φ 2 + 1 + ω x + ω y α k+3 ( 20
)
with initial conditions: α 0 given above,

α 1 = α 0 p φ 2 + 1 + ω x + ω y , α 2 = α 0 2 p φ 2 + 1 + ω x + ω y 2 + p 2 φ 2 2 + 1 + 2pφω x , α 3 = α 0 6 p φ 2 + 1 + ω x + ω y 3 + 3 p φ 2 + 1 + ω x + ω y p 2 φ 2 2 + 1 + 2pφω x + 2 p 3 φ 3 2 + 1 + 3p 2 φ 2 ω x . ( 21
)
The final formula for the convergent series expansion of f and the positivity of the coefficients of this series is stated in the following theorem.

Theorem III.2. The convergent series expansion of f is:

f (ξ) = +∞ k=0 β k ξ k , ( 22
)
where

β 0 = 0, β k+1 = α k (k + 1)! , ( 23)
(α k ) k≥0 are defined in equations [START_REF] Widder | Princeton mathematical series[END_REF], [START_REF] Flajolet | Analytic Combinatorics[END_REF] and β k > 0, for all k ≥ 0.

Proof. The equation ( 23) can be seen as a term-by-term Inverse Laplace Transform, since

L -1 α k λ k+2 = α k ξ k+1 (k + 1)! .
The proper justification for this is based on [START_REF] Widder | An introduction to transform theory[END_REF]Chap. 9]. Moreover, in can be proven using complex analysis arguments [20, Chap. 2.14] that f is an entire function, which implies the convergence of the series in equation [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF]. The positivity of the coefficients β k is deduced from the positivity of the coefficients α k . Indeed, since ϕ(λ) = ∞ k=0 ϕ k λ k is a rational fraction, the general closed-form for its power series coefficients ϕ k can be obtained by hand, based on the partial fraction decomposition of ϕ (or using ratpolytocoeff Gfun procedure):

ϕ k = p k+1 1 + φ k ω x (k + 1) p + φ 2 +          0, k > 0 ω y , k = 0. ( 24
)
It follows that ϕ k > 0, for all k ∈ N. Injecting the series

∞ k=0 ϕ k λ k and ∞ k=0
α k λ k into the differential equation ( 18) and extracting the coefficient of λ k , one gets:

(k + 1)α k+1 = k i=0 ϕ i α k-i , ( 25
)
and hence, since α 0 > 0, it follows that α k > 0, for all k ∈ N and therefore β k > 0, for all k ∈ N.

This entails the final analytic series form of the collision probability formula:

P c = g(R 2 ) = exp -pR 2 +∞ k=0 β k R 2k . ( 26
)
For an efficient computation of truncations of this series, it is important to remark that also the terms β k R 2k satisfy the following linear recurrence:

Proposition III.3. Let c k = β k R 2k . Then P c = exp -pR 2 +∞ k=0 c k , ( 27
)
and (c k ) k≥0 satisfies:

(k + 4)c k+4 = - R 8 p 3 φ 2 ω y (k + 2)(k + 3)(k + 4)(k + 5) c k + R 6 p 2 φ pφ k + 5 2 + 2ω y φ 2 + 1 (k + 3)(k + 4)(k + 5) c k+1 - R 4 p pφ φ 2 + 1 (2k + 5) + φ 2ω y + 3p 2 + ω x + ω y (k + 4)(k + 5) c k+2 + R 2 p(2φ + 1)(k + 3) + p φ 2 + 1 + ω x + ω y k + 5 c k+3 , ( 28
)
with initial conditions

c 0 = α 0 R 2 , c 1 = α 0 R 4 2 p φ 2 + 1 + ω x + ω y , c 2 = α 0 R 6 6 p φ 2 + 1 + ω x + ω y 2 + p 2 φ 2 2 + 1 + 2pφω x , c 3 = α 0 R 8 24 p φ 2 + 1 + ω x + ω y 3 + 3 p φ 2 + 1 + ω x + ω y p 2 φ 2 2 + 1 + 2pφω x + 2 p 3 φ 3 2 + 1 + 3p 2 φ 2 ω x .
Proof. The recurrence for c k is a direct computation (by hand or with Gfun), obtained from the recurrence on α k (see equation ( 20)).

Truncating the series in (27) gives a straightforward algorithm (see Algorithm 1) for computing the probability of collision P c . Besides the fact that this series is the exact analytic formula for the probability, this algorithm features numerical stability. The coefficients c k are positive, so partial sums of the series already provide a lower bound for the exact value of P c . Note that this algorithm requires two exponential function evaluations and O(N ) basic arithmetic operations, where N is the number of terms in the partial sum.

Algorithm 1 Computation of the Probability of Collision.

Input: Parameters: σx, σy, xm, ym; combined object radius: R; Number of terms: N .

Output: Pc -truncated series approximation of Pc.

1: p = 1 2σ 2 y ; φ = 1 - σ 2 y σ 2 x ; ωx = x 2 m 4σ 4
x

; ωy = y 2 m 4σ 4
y ; α0 = 1 2σxσy exp - 1 2 x 2 m σ 2 x + y 2 m σ 2 y ; 2: c0 = α0R 2 ; 3: c1 = α 0 R 4 2 p φ 2 + 1 + ωx + ωy ; 4: c2 = α 0 R 6 6 p φ 2 + 1 + ωx + ωy 2 + p 2 φ 2 2 + 1 + 2pφωx ; 5: c3 = α 0 R
c k+4 = - R 8 p 3 φ 2 ωy (k+2)(k+3)(k+4) 2 (k+5) c k + R 6 p 2 φ pφ(k+ 5 2 )+2ωy φ 2 +1 (k+3)(k+4) 2 (k+5) c k+1 - R 4 p pφ φ 2 +1 (2k+5)+φ(2ωy + 3p 2 )+ωx+ωy (k+4) 2 (k+5) c k+2 + R 2 p(2φ+1)(k+3)+p

C. Computing the probability of collision with guaranteed accuracy

While Algorithm 1 provides already a lower bound for P c , it is useful to also have an upper bound for the truncation error entailed by computing the partial sum Pn = exp -pR 2 n-1 k=0 c k .

The choice of the order n of truncation is strongly dependent upon the trade-off between the numerical complexity involved in the computation and the precision of the obtained result. Thus, it is important to have sharp order-dependent error estimates that will help the user to manage this trade-off. Let us first define what may be interpreted as the truncation of order n:

ε n = P c -Pn (R 2 ), (29)
where Pn = exp -pR

2 n-1 k=0 c k = exp(-pR 2 ) n-1 k=0 α k (k + 1)! (R 2 ) k+1
. The objective of the next proposition is to give accurate bounds on ε n . As will be seen in the following, truncation error bounds also allow for an a priori computation of the truncation order, when a pre-specified precision on the collision probability is required. To avoid rather tedious developments, the proofs of the propositions presented in this Section are gathered in Appendix A.

Proposition III.4. The following bounds hold for the truncation error:

l n ≤ P c -Pn (R 2 ) ≤ u n , ( 30
)
where

l n := α 0 exp(-pR 2 )(pR 2 ) n+1 p(n + 1)! (31)
and

u n := α 0 exp p φ 2 + ωx+ωy p R 2 p 1 + φ 2 + ωx+ωy p R 2 n+1 p 1 + φ 2 + ωx+ωy p (n + 1)! , ( 32
)
for all n ≥ 1.

These bounds are used for computing the probability of collision with a guaranteed accuracy.

Let δ be the probability threshold for executing some collision avoidance maneuver, say δ = 0.001 in practice. The test P c ≤ δ needs to be performed safely, which means that at leastlog 10 δ digits need to be guaranteed in the computation of P c . Practically, for δ = 0.001, 3 or 4 guaranteed digits (after the decimal dot, that is in absolute error) should be sufficient. In the following proposition, a sufficient value for n is given such that the required accuracy holds.

Proposition III.5. Let

N 1 = 2 epR 2 1 + φ 2 + ω x + ω y p , ( 33
)
N 2 =     log 2 α 0 exp pR 2 φ 2 + ωx+ωy p δp √ 2πN 1 1 + φ 2 + ωx+ωy p     , ( 34
)
and

n + 1 = max {N 1 , N 2 } . ( 35
)
Then, u n -l n < δ.

In the proof of Proposition III.4 presented in Appendix A, Proposition A.3 provides two majorant/minorant series ε n and ε n for the truncation error ε n . These can be also used to obtain directly an upper/lower bound for the probability P c as stated below.

Proposition III.6. Let

l 0 = α 0 1 -exp(-pR 2 ) p ,
and

u 0 = α 0 exp p φ 2 + ωx+ωy p R 2 -exp(-pR 2 ) p 1 + φ 2 + ωx+ωy p
.

Then:

l 0 ≤ P c ≤ u 0 . ( 36
)
The formulas obtained for l 0 and u 0 are simple and sufficient in most practical cases. They provide a guaranteed enclosure for the actual value of the probability (see Section V, Table 2 for numerical examples).

Propositions III.6, III.4 and III.5 allow for computing an enclosure of the exact value of the probability with a guaranteed absolute error less then a given threshold δ. This is implemented in Algorithm 2. The advantage of this algorithm is to allow for an a priori computation of a sufficient number n of terms to be taken in the series such that the accuracy requirement holds. This algorithm calls Algorithm 1 (line 10, Algorithm 2) to obtain the sum of the truncated series of n terms i.e., an approximation Pc of the probability of collision. The variables l n and u n (line 11, 12, Algorithm 2)

provide lower and upper bounds on the truncation error. Whence, one obtains an effective enclosure of the true value of the probability in lines 13, 14 of Algorithm 2. Numerical results are given in Section V (see Figure 3 in particular).

Algorithm 2 Computation of probability of collision with guaranteed accuracy.

Input: Parameters: σx, σy, xm, ym; combined object radius: R; Threshold δ.

Output: P c , Pc such that P c ≤ Pc ≤ Pc and Pc -P c ≤ δ.

1: p = 1 2σ 2 y ; φ = 1 - σ 2 y σ 2 x ; ωx = x 2 m 4σ 4
x (n + 1)! ;

; ωy = y 2 m 4σ 4
y ; α0 = 1 2σxσy exp - 1 2 x 2 m σ 2 x + y 2 m σ 2 y ; 2: l0 = α0 1 -exp(-pR 2 ) p ; 3: u0 = α0 exp p φ 2 + ωx+ωy p R 2 -
13: P c = Pc + ln;

14: Pc = Pc + un;

15: return P c , Pc .

16: end if

IV. Comparison with Chan series

In the isotropic case (i.e., σ x = σ y := σ), Chan [5,Chap. 4] reformulates the initial 2-D integral (4) as a one dimension integral over the Rice probability distribution. This integral is then computed analytically in the form of a power series using some knowledge (identities, integral tables) on the modified Bessel function of the first kind and the Rice distribution. We recover the same series applying our procedure described in Section III B. In the isotropic case, the Laplace transform in equation ( 16) is:

L giso (λ) = exp x 2 m +y 2 m σ 2 (2σ 2 λ+1) λ (2σ 2 λ + 1) . ( 37
)
Applying a similar preconditionning method [START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[END_REF][START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], a preconditioning parameter p = 1 2σ 2 is derived. A result similar to Theorem III.2 is easily derived. Note that since σ x = σ y , one has φ = 0,

ω x + ω y = x 2 m + y 2 m 4σ 4 and α 0 = 1 2σ 2 exp - x 2 m + y 2 m 2σ 2 . Define Lfiso (λ) := λ -2 L fiso (λ -1 ), then d Lfiso (λ) dλ = ϕ iso (λ) Lfiso (λ), Lfiso (0) = α 0 , (38)
where

ϕ iso (λ) = ω x + ω y - p pλ -1 . ( 39
)
Lfiso (λ) can be expanded in

+∞ k=0 α k λ k , for all {λ, |λ| > p} and (α k ) k≥0 satisfies (k + 2)α k+2 = -p (ω x + ω y ) α k + (p(k + 2) + ω x + ω y ) α k+1 (40)
with initial conditions: α 0 given above and α 1 = α 0 (p + ω x + ω y ).

Proposition IV.1. The convergent series expansion of the function f iso is:

f iso (z) = +∞ k=0 β k z k (41)
where

β 0 = 0, β k+1 = α k (k + 1)! for all k ≥ 0, (42)
(α k ) k≥0 are defined in equation (40) and β k > 0, for all k ≥ 0.

The advantage in this simpler case is that in (40), one has:

(k + 2) (α k+2 -pα k+1 ) = (ω x + ω y ) (α k+1 -pα k ) , (43)
which easily gives the closed-form for α k+1 -pα k = α0(ωx+ωy) k+1 (k+1)!

, and hence,

α k = α 0 p k k i=0 (ω x + ω y ) i p i i! . ( 44
)
It follows that the final probability becomes:

exp -pR 2 P iso (R 2 ) = ∞ k=0 α k R 2k+2 (k + 1)! = exp - ω x + ω y p ∞ k=0 (pR 2 ) k+1 (k + 1)! k i=0 (ω x + ω y ) i p i i! = exp - ω x + ω y p ∞ i=0 (ω x + ω y ) i p i i! ∞ k=i (pR 2 ) k+1 (k + 1)! ,
or equivalently,

P iso (R 2 ) = exp - ω x + ω y p ∞ i=0 (ω x + ω y ) i p i i! 1 -exp -pR 2 i k=0 (pR 2 ) k k! . ( 45
)
Formula (45) and equivalent ones obtained by simple re-summation are the ones used by Chan, see e.g., [START_REF] Chan | Spacecraft Collision Probability[END_REF]eq. 4.17,Chap. 4], where the notation is

u = R 2 σ 2 and v = x 2 m +y 2 m σ 2 .
In the non-isotropic case, Chan builds an approximation of the integral which allows for the following two changes of variables:

u = r 2 σ 2 x σ 2 y and v = x 2 m σ 2 y + y 2 m σ 2 y in (45)
, used to compute the final approximated probability. In Chan's approach, a closed formula, namely equation (44), is maintained for the sequence α k in the final formula for the non-isotropic case which remains an approximation of the genuine probability of collision. In contrast, the analytic formula obtained in this paper is exact and the sequence α k (eq. ( 20)) is given by a linear recurrence which is numerically stable for a machine computation.

V. Numerical tests

The performance of our method is assessed from two perspectives. First, since it is based on a series expansion, the numerical accuracy varies in function of the number of terms computed.

Algorithm 2 offers an automatic way of computing the number of terms needed for a user-required accuracy. We exemplify it in what follows on practical cases. Second, our method is compared with other methods from the literature concerning the quality of the results obtained. Three algorithms from the literature, namely Alfano's [4], Patera's [START_REF] Patera | General method for calculating satellite conjunction probability[END_REF] and respectively Chan's [START_REF] Chan | Spacecraft Collision Probability[END_REF] have been implemented. The chosen test cases are described in Table 1: the first 12 cases can be found in [5,

Chapter 5] and are supposed to be representative of real short-term encounters; the next 3 cases are real-case scenarios: the data were retrieved from CSMs (Conjunction Summary Messages) sent by the Joint Space Operations Center to the industrial partner of this study; the last two cases were obtained using the physical parameters of test cases number 3 and 5 provided by Alfano [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF]. If more accuracy is needed, the subsequent lines of Algorithm 2 are employed. Figure 3 shows the number n of series terms needed in Algorithm 2 for a requested accuracy δ ranging from 10 -1 to 10 -13 . We observe that the number of terms for the absolute error to reach the machine precision (10 -13 ) is less 40 in all of Chan's Cases and CSM cases. Note that this number is computed a priori and it is a sufficient number. The actual number of terms for the accuracy to be met may be smaller in some cases, but for these practical examples it is not conservative. Chan's cases 8 and 10 are not drawn since zero terms (that is the values computed in Table 2 with the simple bounds) also performs well, but it fails for very low probabilities like in test cases 8 and 10. On the other hand, Chan's method gives non negligible relative errors for some cases -namely 1, 2, 3, 4, 8, 10, 12 and CSM1-3. As far as precision is concerned, it is definitely the least effective. It is not surprising since it is based on an additional approximation with respect to the original short-term encounter model. For that reason, it gives meaningless results for the two test cases provided by Alfano [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF], see Table 5. These examples were originally designed to compare the efficiency of several methods of the literature and are somehow more tedious as far as computation is concerned. They are challenging also for our method, since the number of terms to be considered in the series expansion is important.

Nevertheless, the new method gives satisfactory results.

Concerning timings, our method is very fast: for each case tested, the results are obtained in less than one second; in frequent cases, when the bounds l 0 and u 0 are sufficient (lines 2 -3, Algorithm 2), the response is almost instantaneous (10 -5 seconds).

Case Collision Probability (-)

VI. Conclusions

A new method to compute the probability of collision for short-term space encounters between two spherical objects under Gaussian-distributed uncertainty has been proposed. The formula has the form of a product between an exponential term and a convergent power series with positive coefficients. It is derived by use of Laplace transform theory and D-finite functions properties. The variable in the power series is the squared radius of the combined object. The series coefficients depend on the other parameters of the encounter, namely the miss distance and the standard deviations of the relative position of the objects in the encounter plane. It was shown that the formula is in fact an improvement on the only other analytical method of the literature -namely Chan's -since it is more general. Compared with state-of-the-art numerical methods, this work has several advantages. First of all, it provides an accurate algorithm to compute the probability of collision based on an analytic expression. Moreover, it gives a fast way to compute an estimate of the risk, which is often close-enough to the real value for the decision maker.

As a short term objective, we intend to address, in a finer manner, the cases for which the number of terms to be considered in the series expansion is too high. We will use for that purpose, a saddle-point based method [21, Chapter VIII]. In addition, a parameter sensitivity analysis similar to the one done by Chan, would be useful since we expect to obtain similar results (although with a better accuracy on the probability, due to the advantages of our method). In order to prove the Proposition III.4, some preliminary results are needed. First, remember that

ε n = P c -Pn (R 2 ) = P c -exp(-pR 2 ) n-1 k=0 β k (R 2 ) k+1 = P c -exp(-pR 2 ) n-1 k=0 α k (k + 1)! (R 2 ) k+1 . (A1)
Then, recall that (α k ) k≥0 is known only through its recurrence and a closed formula is unlikely to be found. So, in order to compute effective truncation error bounds, two sequences (α k ) k≥0 and (α k ) k≥0 are used, which have closed-forms and provide respectively upper and lower bounds for Proof. The proof follows by induction. First, one has α 1 = α 0 ϕ 0 = α 0 p 1 + φ 2 + ωx+ωy p , so

α 1 = α 0 p ≤ α 0 ϕ 0 ≤ α 0 p 1 + φ 2 + ω x + ω y p = α 1 . (A2)
To prove the induction step, let us first note:

α 0 ϕ k α k+1 = 1 + φ k ω x (k + 1) p + φ 2 +          0, k > 0 ω y p , k = 0. ( A3
)
and

α 0 ϕ k α k+1 =                1 + φ k ω x (k + 1) p + φ 2 1 + φ 2 + ωx+ωy p k+1 , k > 0 1, k = 0, (A4)
which implies:

α 0 ϕ k α k+1 ≥ 1 ≥ α 0 ϕ k α k+1 . (A5)
To prove these two inequalities, remember that φ = 1 -σ 2 y σ 2

x ≥ 0 since σ 2 x and σ 2 y are assumed to be respectively the largest and smallest eigenvalues of Σ xỹ . Indeed, α 0 ϕ k α k+1 ≥ 1 is easily obtained from the positivity of all the terms (φ, φ k , ω x , ω y , p) involved in (A3). The second inequality comes from the following series of inequalities and from the binomial formula, still using positivity of the previous terms:

ε n (R 2 ) ≤ P c -Pn (R 2 ) ≤ ε n (R 2 ), ( A9
)
with the particular case P0 (t) := 0, for n = 0.

Proof. From equation (26) and β 0 = 0, one gets: P c -Pn (R 2 ) = exp(-pR 2 ) ∞ k=n α k (k + 1)! R 2(k+1) .

Apply then Corollary A.2.

To conclude the proof of the Proposition III.4, note that in general, for γ > 0,

1 γ (γt) n+1 (n + 1)! ≤ ∞ i=n γ i t i+1 (i + 1)! ≤ e γt γ (γt) n+1 (n + 1)! , ( A10
)
which gives straightforwardly ε n (R 2 ) ≥ l n and ε n (R 2 ) ≤ u n , hence (30).

Proof of Proposition III.5

From equations (31), (32) and using the Stirling inequality [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF] on l n , u n -l n = l n exp pR 2 .

For n + 1 = N 1 in the previous inequality, the obtained bound for the absolute error is:

u N1-1 -l N1-1 < 1 2 N1 α 0 exp pR 2 φ 2 + ωx+ωy p p √ 2πN 1 1 + φ 2 + ωx+ωy p . ( A11
)
If this is not less than δ, one has for all N ≥ N 1 :

u N -1 -l N -1 ≤ 1

Figure 1 :

 1 Figure 1: Encounter plane and frame (e x, e ỹ , e z )

Figure 2 :

 2 Figure 2: Rotation to the principal axis of the covariance matrix in the encounter plane

10 :

 10 for k ← 0 to N -1 do 11: s ← s + c k ;12: end for13: Pc ← exp -pR 2 s; 14: return Pc.

3 Chan 2 9 . 3 Chan 3 6 . 3 Chan 4 6 .

 393636 742 × 10 -3 9.741 × 10 -3 9.754 × 10 -3 9.742 × 10 -3 9.742 × 10 -181 × 10 -3 9.181 × 10 -3 9.189 × 10 -3 9.181 × 10 -3 9.181 × 10 -571 × 10 -3 6.571 × 10 -3 6.586 × 10 -3 6.571 × 10 -3 6.571 × 10 -125 × 10 -3 6.125 × 10 -3 6.135 × 10 -3 6.125 × 10 -3 6.125 × 10 -3Chan 5 1.577 × 10 -5 1.577 × 10 -5 1.577 × 10 -5 1.577 × 10 -5 1.577 × 10 -5 Chan 6 1.011 × 10 -5 1.011 × 10 -5 1.011 × 10 -5 1.011 × 10 -5 1.011 × 10 -5

(α k ) k≥0 . Proposition A. 1 .

 1 Let α k = α 0 p k and α k = α 0 p k 1 + φ 2 + ωx+ωy p k . Then α k ≤ α k ≤ α k , for all k ∈ N.

  exp(-pR 2 )

			p 1 + φ 2 +	p ωx+ωy	;
			ωx+ωy p	;
	8: N2 = log 2	α0 exp pR 2 φ 2 + δp √ 2πN1 1 + φ 2 + ωx+ωy p p ωx+ωy	;
	11: ln =	α0 exp(-pR 2 )(pR 2 ) n+1 p(n + 1)!	;
	12: un =	α0 exp p φ 2 + p 1 + φ ωx+ωy p 2 + R 2 p 1 + φ 2 + ωx+ωy p	ωx+ωy p	R 2 n+1

4: if u0 -l0 ≤ δ then 5: return [l0, u0]; 6: else 7: N1 = 2 epR 2 1 + φ 2 + 9: n = max {N1, N2} -1 10: Pc ← Algorithm 1 (σx, σy, xm, ym, R, n);

Table 1 :

 1 Inputs for test cases from Chan[START_REF] Alfriend | Probability of collision error analysis[END_REF][START_REF] Akella | Probability of collision between space objects[END_REF][START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Coppola | Including Velocity Uncertainty in the Probability of Collision Between Space Objects[END_REF][START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF][START_REF] Foster | A parametric analysis of orbital debris collision probability and maneuver rate for space debris[END_REF][START_REF] Alfano | Aerospace support to space situation awareness[END_REF][START_REF] Chan | Collision probability analyses for earth-orbiting satellites[END_REF][START_REF] Lasserre | Solving a class of multivariate integration problems via laplace techniques[END_REF][START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF], CSMs (1-3) and Alfano (3,5) Algorithm 2) provide enough accuracy in most practical cases, as shown in Table2. Except in Alfano's cases, they seem to be very effective due to their simplicity and may suffice in most practical cases.

	Case			Input parameters (m)
	#	σx	σy	R	xm	ym
	Chan 1	50	25	5	10	0
	Chan 2	50	25	5	0	10
	Chan 3	75	25	5	10	0
	Chan 4	75	25	5	0	10
	Chan 5	3,000	1,000	10	1,000	0
	Chan 6	3,000	1,000	10	0	1,000
	Chan 7	3,000	1,000	10	10,000	0
	Chan 8	3,000	1,000	10	0	10,000
	Chan 9	10,000	1,000	10	10,000	0
	Chan 10	10,000	1,000	10	0	10,000
	Chan 11	3,000	1,000	50	5,000	0
	Chan 12	3,000	1,000	50	0	5,000
	CSM 1 152.8814468961533 57.918666623295984 10.3 60.583685340533115 84.875546447209487
	CSM 1 5,756.840725983703 15.988242371297744 1.3 115.0558998093139 -81.618369910317043
	CSM 3 643.4092722122279 94.230921098486149 5.3 693.4058939950484 102.1772470067133
	Alfano 3 114.2585190378857 1.410183033040157 15 0.159164620813659 -3.887207383647396
	Alfano 5 177.8109003935867 0.037327944173609 10 2.123006718041866 -1.221789517557463
	A. Analysis of our method
	In Proposition III.6 two very simple formulas (evaluations of exponentials) are provided for
	guaranteed lower and upper bounds of the probability of collision. These two simple formulas (line
	2, 3,

Table 2 :

 2 Guaranteed correct digits obtained with formulas in Proposition III.6. The notation e.g. 0.15[START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF] 65] • 10 -8 stands for the fact that the probability is surely enclosed in the interval [0.1518 • 10 -5 , 0.1565 • 10 -8 ] and helps visualizing the guaranteed digits.

  Chan 7 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -8 10 -27 3.216 × 10 -27 3.219 × 10 -27 3.219 × 10 -27 Chan 9 3.033 × 10 -6 3.033 × 10 -6 3.033 × 10 -6 3.033 × 10 -6 3.033 × 10 -6 10 -28 9.645 × 10 -28 9.656 × 10 -28 9.656 × 10 -28 Chan 11 1.039 × 10 -4 1.039 × 10 -4 1.039 × 10 -4 1.039 × 10 -4 1.039 × 10 -4Chan 12 1.564 × 10 -9 1.564 × 10 -9 1.556 × 10 -9 1.564 × 10 -9 1.564 × 10 -9

	Chan 8 3.219 × Chan 10 0 0 9.656 ×

Table 3 :

 3 Comparison of collision probability value -with 4 significant digits-for Chan's test cases number 1 to 12. The digits different from the reference value are represented in bold. .9002 × 10 -3 1.9001 × 10 -3 1.8934 × 10 -3 1.9002 × 10 -3 1.9002 × 10 -3 CSM 2 2.0553 × 10 -11 2.0552 × 10 -11 2.0135 × 10 -11 2.0553 × 10 -11 2.0553 × 10 -11 CSM 3 7.2004 × 10 -5 7.2000 × 10 -5 7.2000 × 10 -5 7.2003 × 10 -5 7.2003 × 10 -5

	Case		Collision Probability (-)
	Alfano	Patera	Chan	Algorithm 2	Reference
	CSM 1 1

Table 4 :

 4 Comparison of collision probability value -with 5 significant digits-for tests cases from CSMs. The digits different from the reference value are represented in bold. No. 3 1.0038 × 10 -1 1.0087 × 10 -1 3.1264 × 10 -2 1.0038 × 10 -1 1.0085 × 10 -1 Alfano's No. 5 4.4712 × 10 -2 4.4520 × 10 -2 1.6618 × 10 -77 4.4509 × 10 -2 4.4499 × 10 -2

	Case		Collision Probability (-)
	Alfano	Patera	Chan	Algorithm 2	Reference
	Alfano's

Table 5 :

 5 Comparison of collision probability value -with 5 significant digits -for test cases from[START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF].

  It is then quite immediate to get the following result.Proof. Direct application of Proposition A.1 and Theorem III.2.Upper and lower bounds for the truncation errors are now obtained based on the above minorant/majorant sequences.

	Proposition A.3. Let Pn (t) = exp(-pt) ε n (t) = e -pt ∞ i=n α i (i + 1)! t	n-1 i=0		α i (i + 1)!	t i+1 , ε n (t) = exp(-pt)	∞ i=n	α i (i + 1)!	t i+1 and
		1 + φ k ω x (k + 1) p	+	φ 2			1 +	ω x (k + 1) p	+	φ 2
										≤ 1 +	(ω x + ω y )(k + 1) p	+	φ 2
										≤ 1 +	(ω x + ω y ) p	+	φ 2
											n
											ϕ k α n-k	(A6)
											k=0
	Since α n-k =	α 0 α n+1 α k+1	and α n-k =	α 0 α n+1 α k+1	, it comes that:
		α n+1	n k=0	α 0 ϕ k α k+1	≤		n k=0	ϕ k α n-k ≤ α n+1	n k=0	α 0 ϕ k α k+1	.	(A7)
	Hence,
								n
			α n+1 (n + 1) ≤		ϕ i α n-i ≤ α n+1 (n + 1),	(A8)
								i=0
	and the conclusion follows from (25).
	Corollary A.2. The coefficients (β k ) k≥0 satisfy	α k (k + 1)!	≤ β k+1 ≤	α k (k + 1)!	.

k+1

.

Suppose α k ≤ α k ≤ α k , ∀k ≤ n. Then: n k=0 ϕ k α n-k ≤ n k=0 ϕ k α n-k ≤ i+1 . For each n ≥ 0, one has

  [START_REF] Akella | Probability of collision between space objects[END_REF]

		N α 0 exp pR 2 φ 2 + p √ 2πN 1 1 + φ 2 + ωx+ωy p p ωx+ωy	.	(A12)
	Finally, one has (by construction):
	1 2	N2 α 0 exp pR 2 φ 2 + p √ 2πN 1 1 + φ 2 + ωx+ωy p p ωx+ωy	≤ δ,		(A13)

so, n + 1 = max {N 1 , N 2 } is suitable.

Acknowledgments The Authors would like to thank Alexandre Falcoz from Airbus Defence and Space and the ANR FastRelax Project for the grants that partly supports this activity.

are sufficient for the whole range of absolute error considered.

On the other hand, for Alfano's cases, our Algorithm needs much more terms. One observes in Table 2 that the bounds obtained for these two examples are meaningless. In order to reach a good accuracy, Algorithm 2 computes n = 689 for Alfano's Case 3 and n = 10 13 for Alfano's Case 5. This is conservative, since 90 terms for Case 3 and respectively 37000 terms for Case 5 are sufficient to obtain the value given in Table 5. This shows however that in some degenerate cases the number of terms needed may increase drastically. 1.

B. Comparison with other methods

The corresponding results for the probability of collision obtained with different methods are summarized in Tables 3, 4 and5. The reference values in Table 3 were provided by NASA [5] using Foster's method. For Table 4, they were given by the industrial partner and for Table 5 In general, for γ > 0,

and:

By multiplying the above expressions by e -pR 2 , one obtains: l 0 = ε 0 (R 2 ) = α 0 exp(pR 2 ) -1 p and:

Similarly, one gets: ε 0 (R 2 ) + Pn (R 2 ) ≤ ε 0 (R 2 ).


Clement Poncelet

Florent Jacquemard

Offline methods of conformance testing for Antescofo

Keywords: Model-based testing, offline, Antescofo tests, manual Méthodes test basé sur model, différé, tester Antescofo

This report presents procedures for testing Antescofo, a score-based Interactive Music System (IMS). It describes common Model-Based Test (MBT) vocabulary and presents methodology and use cases.

Model Based Testing

Figure 1 depicts in its upper half a reactive system's Implementation Under Test (IUT) interacting with an environment RealENV, and in its lower half, two formal specifications of the latter, resp. S and E.

The environment model E is a formal description of what can be expected from the environment. In our case, it is the definition of the set of all possible input traces, i.e. all the potential interpretations of musicians to be tested on a given mixed-score.

Note that since IMS are realtime systems, we need to express time in E and S, like in: "one message m has to be emitted one beat after the first event e 1 of the musician".

We present the application of model based testing (MBT) techniques to a score based IMS called Antescofo. Roughly, our method proceeds with the following steps, depicted in Figure 2. First, a given mixed score is compiled into an intermediate representation (IR). This formalism can be presented as a table of finite state machines extended with delays and asynchronous communications with the environment (using input and output symbols) and between machines. It is viewed as the model M describing Antescofo reactions according to the score. Note that the model contains the environment (E) and the system (S) specifications.

The IR is then transformed into a timed automata (TA) network, assuming some restrictions. TA's are processed by tools from the Uppaal suite to generate test cases or expected output suites. In particular we use CoVer in some scenarios to make covering suites of test cases: input timed traces together with the corresponding expected output.

The input traces are then sent to the IUT (as a virtual musician's performance of the score) and the real outcome of the IUT is analyzed, resulting in a test verdict. Figure 3 describes roughly the architecture of Antescofo, which is made of two main modules. A listening machine (LM) decodes an audio or midi stream incoming from a musician and infers in realtime: (i) the musician's position in the given mixed score, (ii) the musician's instantaneous pace (tempo, in beats per minute). These values are sent to a reactive engine (RE) which schedules the electronic actions to be played, as specified in the mixed score. The actions are messages emitted on time to an audio environment. Therefore, the information exchanged between LM and RE as well as between RE and the output environment of the system is made of discrete events.

IO and Design Specific Language representation

Let O be a set of output messages (also called action symbols and denoted a) which can be emitted by the system and let I be a set of event symbols (denoted e) to be detected by the LM (i.e. positions in score).

In our Antescofo's DSL representation, an action is a term act(d, s, al ) where d is the delay before starting the action, s is either an atom in O or a finite sequence of actions (such a sequence is called a group), and al is a list of attributes.

A mixed score is a finite sequence of input events of the form evt(e, d, s) where e ∈ I, d is a duration and s is the top-level group triggered by e. Sequences We consider here two time units for expressing delays and durations d: (i) the number of beats (default unit): a logical time unit traditionally used in music scores that we call relative time, and (ii) milliseconds (ms), referred to as physical time.

The reconciliation of the relative and physical times is done through the detected tempo values. We consider a black-box testing conformance approach for Antescofo, based on timed traces comparisons and assume a given mixed score with a default tempo value.

Timed traces

A timed trace is a sequence of pairs a, t made of a symbol a ∈ I ∪ O, and a timestamp t ∈ R + , either in physical or relative time.

A trace containing symbols exclusively in I (resp. O) is called an input trace (resp. an output trace). We denote below T in (resp. T out ) the set of input (resp. output) traces with relative timestamps.

The ideal trace is the input trace consisting in the projection of all events in the mixed-score with their duration.

By definition of music performance, traces of real executions can be arbitrarily far from ideal traces: the tempo and delays can diverge from the written values (the musician adding her/his own expressiveness values), and moreover the musician can do errors during a performance (missing notes).

A test case is a pair t in , t out ∈ T in × T out where t in ∈ E and t out = S(t in ). t out = S(t in ) is a simulation of t in on the specification returning t out . We note too t out = IUT(t in ) an Antescofo execution with t in as input trace and t out as monitored real trace. We recall that: The principle (Figure 5) is to execute a real implementation under test (IUT) in a testing framework. When the source code of the IUT is not known and only its input and output are observed, we call it black-box testing. In conformance model-based testing (MBT), a formal specification, or model M of the system is required to generate automatically some test data. This comprises input test data t in , sent to the IUT, and theoretically expected output test data t out , computed from t in using S. The latter t out is then compared to the real output test data t out , obtained from the IUT when it receives t in , in order to produce a test verdict.

Execution

The execution of a test case t in , t out consists in several tasks summarized in the following definition of conformance. First, we pass the events of t in to the IMS Antescofo, respecting the timestamps. Second, we monitor the outcome of the IUT in an output trace t out = IUT(t in ). Finally, we compare t out to t out = S(t in ). We define the conformance of the IUT to S wrt E as:

∀t in ∈ E, IUT(t in ) = S(t in ).
This is a particular case of the relation rtioco e .

This definition makes sense only if the timestamps of t in , t out and t out are in the same time unit. We will show how this important issue is addressed in practice with several options for the conversion of all traces into physical time (thanks to the addition of tempo values).

One trace format is defining to link test softwares and Antescofo symbols. For input traces in use cases, we use the i, d, t input form (relations between a, t and i, d, T are explained section 4).

Tests use cases

Next sub sections present different use cases, implemented to test Antescofo. Each case has a proper handling of the MBT method which can be divided in three independent parts:

• The input, the execution and the conformance.

Input: From the inputs choice, the tester can test her/his mixed-score on a large scale or an exhaustive set of input traces, or just with a specific performance. The first case uses the automatic input trace generation with Uppaal and CoVer or the fuzz method (depending of the length of the mixed-score), and the latter case needs a single input trace or an audio file performance.

Execution: The execution scenario chosen, for a test, defines the black box bounds. Four scenarios are implemented testing from the RE only to the entire IUT in its MAX patch version.

Conformance: The conformance compares two traces. One (the compared trace) has to be an Antescofo output with physical timestamps and their corresponding tempos. So the second (expected trace) can be in relative or physical time. We can choose here whether the real and expected traces are conform or several real traces (obtained from different mixed scores) provoke the same Antescofo behavior (i.e all are conformed to one excepted trace).

To continue the three first use cases concern different scenarios of execution. The next presents the fuzz solution and the multiple conformances case is shown last. The messages sent by RE are traced in t out , with timestamps in physical time (this functionality is built in the current RE). Finally, the timestamps in t out are converted from relative time to physical time using the tempo values in t in , in order to be comparable to t out . In this scenario, the IUT is the RE (the LM is idle).

Testing the RE

Testing the RE with tempo detection

In this second scenario the tempo values T are not read in t in by the adapter but instead inferred by the LM (the adapter is calling an appropriate method in LM). The rest of the scenario is similar to Section 3.2.1. The values of detected tempo are stored by the adapter and used later to convert the dates in the expected trace t out from relative to physical time. In this case, the IUT is somehow the RE plus the part of the LM in charge of tempo inference.

Testing the whole system

This scenario is the most general. It is executed in a version of Antescofo embedded into MAX (as a MAX patch), using an adapter which is another MAX patch. The adapter iteratively reads triples i, d, T in a file containing t in , and plays them as MIDI events, using the duration d converted to physical time with T . The audio stream generated is then sent to the LM, and the output of the RE is traced in t out as before.

Note that here, the RE uses the tempo values detected by the LM, which will differ from the tempo values in t in . Therefore, the former are saved by the adapter (in MAX the detected tempo is available as an outlet of the antescofop atch), and used later to convert the dates in t out from relative to physical time. In this realistic scenario, the IUT is the whole Antescofo system.

In an alternative scenario, the adapter ignores the tempo values in t in and uses tempo detected by the LM, similarly to Section 3.2.2. Note that in these scenarios, the tests are executed in real-time and not in a fast-forward mode like in Sections 3.2.1.

Fuzz testing

The problem handled here is the scalability. As explained before Uppaal tools as CoVer, used to generate a covering suite of input traces, can be overwhelming when too long mixed-scores are given in Another generation can be done by fuzzing the ideal trace and set fuzzed performances in the environment (E) of the model. The output after simulating E on the specification S of the mixed-score is the corresponding Antescofo reaction (the expected trace t out ).

Fuzzing can be performed by an external program which randomly chooses (in a given rate) a different value to the ideal trace. Note that inputs can be manually given by testers and/or composers (i.e via an Ascograph 1 extension). It will not change the rest of the test scenario since the performance is translated in Uppaal environment automaton E.

Testing from an audio file

The problem handled here is the model sub-language restriction. Our model doesn't cover the entire possibilities of the Antescofo DSL. However the verdict is done on the output traces of Antescofo.

This scenario needs a specifiable mixed-score named reference, an audio recorded file and a mixed-score folder. The idea is to derive once or more times a complex mixed-score (which can't be specified) and give a simpler version (reference) of the same score that can be specified (for a performance or a set). In a way it is an Antescofo specification.

1 Ascograph is an Antescofo tool software to manage Antescofo and create a mixed-score graphically, http://repmus.ircam.fr/mutant/ascograph. The audio represents the input trace which will be converted into t in and translated into E. The reference is translated into the specification (S) to generate the expected trace with E. From here, the scenario takes the expected trace and for each mixed-score in the folder (named Tested Mixed-scores) checks the conformance of the result of Antescofo with the audio input.

Note that folder's mixed scores have no more DSL restriction and have to trace the same Antescofo behavior. The conversion of audio into t in is done by passing the audio file to Antescofo which creates the corresponding input trace (that is detected by Antescofo).

Trace format and example

An example is depicted with the revisited Einspielung mixed-score. Only the event and action labels are added of the Nune's original version (we need single labels to identify throw severals programs each event and action of the tested mixed-score). An extract from Ascograph test-editor is shown Figure 11.

Input traces

An input trace is a list of input symbols, i.e mixed-score events. The format used during tests is a suite of lines representing the played events and containing:

• the identifier (notenum)

• the duration of the event (in beat)

• the current tempo during the duration (in BPM)

• the mixed-score label

• the list of notes playing for this event

The two latest fields are optional, and the tempo can be omitted for a relative input trace. As an example, Figure 12 depicts an absolute input trace, i.e an Einspielung performance.

In absolute cases, we represent the theoretic a, t with i, d, T , i is the event identifier of Antescofo and d, T the duration and tempo used to compute the time stamp t in absolute time (t = t -1 + d with d = d.6.10 4 T ).

Output traces

Output traces are expected and real traces. To facilitate the understanding and further handles, we put too event detections in our traces (which aren't just output traces). Their general form is a list of time-stamped symbols.

Expected traces

An expected trace is a trace obtained by the Specification after a simulation of an input trace. Since it comes from the model we choose to keep all the timestamps in relative time unit (beats). Events are kept for testing the detection or the difference between a simulation and an Antescofo trace that can be compared latter. The events are listed with:

• a keyword event • its label

• its duration (in beat)

• its time-stamp (in beat)

• its related tempo And actions with:

• its time-stamp (in beat)

• its label All the time-stamps are taken from the begging of the trace (a computation is needed to have the time elapsed since the nearest passed event to translate, with the corresponding tempo, the good absolute duration and time-stamp). An example is shown Figure 13.

Real traces

A real trace is a trace obtained with Antescofo test trace option and after an execution of an input trace for the tested mixed-score. Information is taken from Antescofo variables during the execution. An event has: • its identifier (notenum)

• its label

• its absolute time-stamp (gettime, the current time [but currently from a virtual clock])

• its dynamic score position (rnow)

• its current detected tempo (rt tempo)

A message:

• its label

• its absolute time-stamp (gettime)

• its dynamic score position (rnow)

• the current tempo Be careful the used time-stamps here (from gettime) are in ms. The real trace example is reported Figure 14 5 Verdict For the moment the verdict was only used for the development of the test method but aims at being a complete representation of test information.

The verdict version reports in a file:

• the version and options of the test

• the input trace linked

• a table of real/expected symbols

• a verdict, ok or ko if the test pass or not

In the table is annotated for each symbol the label, its absolute time-stamp and its relative one (in []). The current tempo for events *, missed symbols (x for messages, o for events). And the delta when the difference of the real/expected time-stamps is greater than the test tolerance.

Figure 15 and 16 present two extracts of the verdict of Einspielung (the beggining and the end).

START of traces

EVENT e0 0.28656 (at 0) tempo=64.2 0 a0 0 a1 0 a2 0 a3 0 a4 0 a5 0 a6 0 a7 0 a8 0 a9 0 a10 0 a11 0 a12 0 a13 0 a14 0 a15 0 a16 0 a17 0.14285 a18 0.14285 a19 0.14285 a20 0.14285 a21 0.14285 a22 0.14285 a23 0.2857 a24 0.2857 a25 0.2857 a26 0.2857 a27 0.2857 a28 0.2857 a29 EVENT e1 0.14199 (at 0.28656) tempo=57.5 0.42855 a30 0.42855 a31 0.42855 a32 0.42855 a33 0.42855 a34 0.42855 a35 EVENT e2 0.28999 (at 0.42855) tempo=60 0.5714 a36 0.5714 a37 * T: 92.5 BPM |----------------------------------------------------------------| Error :: Test KO

--> 0.28656 --> e1_1 --> 0.14199 --> e1_2 --> 0.28999 --> e1_3 --> 0.443 --> e1_4 -- __________________________________________________________________ |- Antescofo Trace Expected Trace -| |- TimeStamp [Estimate Beat] TimeStamp [Beat] -| |----------------------------------------------------------------| a1_0_0 0 [ 0] a0 0 [ 0] a1_0_1 0 [ 0] a1 0 [ 0] a1_0_2 0 [ 0] a2 0 [ 0] a1_0_3 0 [ 0] a3 0 [ 0] a1_0_4 0 [ 0] a4 0 [ 0] a1_0_5 0 [ 0] a5 0 [ 0] a1_0_6 0 [ 0] a6 0 [ 0] a1_0_7 0 [ 0] a7 0 [ 0] a1_0_8 0 [ 0] a8 0 [ 0] a1_0_9 0 [ 0] a9 0 [ 0] a1_0_10 0 [ 0] a10 0 [ 0] a1_0_11 0 [ 0] a11 0 [ 0] a1_1_1 0 [ 0] a12 0 [ 0] a1_1_2 0 [ 0] a13 0 [ 0] a1_1_3 0 [ 0] a14 0 [ 0] a1_1_4 0 [ 0] a15 0 [ 0] a1_1_5 0 [ 0] a16 0 [ 0] a1_1_6 0 [ 0] a17 0 [ 0] * Mesure1 0 [ 0] e0 0

  Figure 1: Specification : reality (top) and models (down)

Figure 3 :

 3 Figure 3: Architecture of Antescofo

  Figure 4: Grammar of the handled AST-mixed score

Figure 5 :

 5 Figure 5: Score-based IMS testing workflow with traces

Figure 8 :

 8 Figure 8: Testing scenarios of Section 3.2.3.

  Figure 10: audio scenario case

Figure 11 :

 11 Figure 11: Einspielung test version extract in the Ascograph text-editor

Figure 12 :

 12 Figure 12: An Einspielung input trace extract

Figure 13 :

 13 Figure 13: An Einspielung expected trace extract

Figure 16 :

 16 Figure 16: An Einspielung verdict extract (last lines)

RR n°8700

RR n°8700 RESEARCH CENTRE PARIS -ROCQUENCOURT Domaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399