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 Power Series Expansion based 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 . 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.

For manipulating D-finite functions/P-recursive sequences, we use the Maple Gfun (version 3.65) package [START_REF] Salvy | Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable[END_REF]. This allows us to show that g is D-finite and to obtain the recurrence formula satisfied by the coefficients of its power series * . This recurrence formula can of course be also obtained/checked in a rather tedious way by pen and paper derivation.

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. * For convenience, the Maple files used are available at http://homepages.laas.fr/mmjoldes/CollisionProba/

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 which are detailed in Appendix B). 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

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] Salvy | Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable[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 [20, Chap. 9], [21, 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 | An introduction to transform theory[END_REF], [START_REF] Widder | Princeton mathematical series[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 [21, Chap. 2.14] that f is an entire function, which implies the convergence of the series in equation [START_REF] Flajolet | Analytic Combinatorics[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 [START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[END_REF] 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 [START_REF] Widder | An introduction to transform theory[END_REF]).

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 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 ; 6: for k ← 0 to N -5 do 7: 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
)
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 ≤ δ. (n + 1)! ;

1: p = 1 2σ 2 y ; φ = 1 - σ 2 y σ 2 x ; ωx = x 2 m 4σ 4
13: P c = Pc + ln;

14: Pc = Pc + un;

15: return P c , Pc .

16: end if 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 , (

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 [START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF], 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) 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 and 5. 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 they were obtained from Monte Carlo trials. All tests were performed with Matlab c R2014a on an Intel Xeon at 3.60GHz.

Since in our method the required accuracy can be set a priori in Algorithm 2, the obtained values are identical (in most cases) or very close to the reference. For Chan's test cases, number 1 to 12, Patera's method gives also 0% of relative error. For the same examples, Alfano's method 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). Table 5: Comparison of collision probability value -with 5 significant digits -for test cases from [START_REF] Alfano | Satellite conjunction monte carlo analysis[END_REF].

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 [START_REF] Flajolet | Analytic Combinatorics[END_REF]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

(α k ) k≥0 . Proposition A.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.
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:

1 + φ k ω x (k + 1) p + φ 2 ≤ + ω x (k + 1) p + φ 2 ≤ 1 + (ω x + ω y )(k + 1) p + φ 2 ≤ 1 + (ω x + ω y ) p + φ 2 k+1 . Suppose α k ≤ α k ≤ α k , ∀k ≤ n. Then: n k=0 ϕ k α n-k ≤ n k=0 ϕ k α n-k ≤ n k=0 ϕ k α n-k (A6) 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+1 (n + 1) ≤ n i=0 ϕ i α n-i ≤ α n+1 (n + 1), (A8)
and the conclusion follows from [START_REF] Rubel | Entire and meromorphic functions[END_REF].

It is then quite immediate to get the following result.

Corollary A.2. The coefficients (β k ) k≥0 satisfy

α k (k + 1)! ≤ β k+1 ≤ α k (k + 1)! .
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-1 i=0 α i (i + 1)! t i+1 , ε n (t) = exp(-pt) ∞ i=n α i (i + 1)! t i+1 and ε n (t) = e -pt ∞ i=n α i (i + 1)! t i+1 . For each n ≥ 0, one has ε 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 1 + φ 2 + ω x + ω y p 1 + φ 2 + ω x + ω y p n -1 < l n exp pR 2 1 + φ 2 + ω x + ω y p 1 + φ 2 + ω x + ω y p n < α 0 exp(-pR 2 )(epR 2 ) n+1 p 2π(n + 1)(n + 1) n+1 exp pR 2 1 + φ 2 + ω x + ω y p 1 + φ 2 + ω x + ω y p n <   epR 2 1 + φ 2 + ωx+ωy p n + 1   n+1 α 0 exp pR 2 φ 2 + ωx+ωy p p 2π(n + 1) 1 + φ 2 + ωx+ωy p .
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 2 N α 0 exp pR 2 φ 2 + ωx+ωy p p √ 2πN 1 1 + φ 2 + ωx+ωy p . (A12)
Finally, one has (by construction): In order to avoid this cancellation problem, for this simple function, one could use the series of exp(x) for which all terms are positive and compute afterwards 1 exp(x) . This is the most basic instance of the idea employed in [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] that consists in evaluating two "well conditioned" series, say f and ψ such that g = f ψ . In this example, ψ = exp and f = 1.

1 2 N2 α 0 exp pR 2 φ 2 + ωx+ωy p p √ 2πN 1 1 + φ 2 + ωx+ωy p ≤ δ, (A13) so, n + 1 = max {N 1 , N 2 } is suitable. Example B.1. A classical example of catastrophic cancellation is for g(ξ) = exp(-ξ). One has g(ξ) = ∞ i=0 (-1) i i! ξ i . The values (-1) i i! ξ i
This idea is developed in [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], by Gawronski, Müller and Reinhard, who provide a method (called GMR from now on) for evaluating an entire function on some complex sector with reduced cancellation. We use a simple instance of their method in the following, since we focus only on positive real line evaluation. Before explaining GMR method and how to apply it for our case, we review some classical complex analysis notions (see [START_REF] Levin | Lectures on Entire Functions[END_REF]Chap.1,[START_REF] Alfano | Aerospace support to space situation awareness[END_REF] for example).

GMR Method in a nutshell

a. Some classical complex analysis notions. Let g(ξ) = ∞ n=0 g n ξ n be an entire function (i.e., analytic on the whole complex plane, lim

n→∞ n |g n | = 0).
• The order (of growth) ρ of g is given by ρ = lim sup The order relates to the coefficients:

ρ = lim sup n→∞ n log n log |g n | -1 . (B3) • If 0 < ρ < ∞, the type σ is given by σ = lim sup n→∞ log M g (r) r ρ . ( B4
)
The following relation between order, type and coefficients growth holds: if g is of order ≤ ρ:

σ = 1 eρ lim sup n→∞ n |g n | ρ/n . (B5)
• The indicator function with respect to ρ is h:

h(θ) = lim sup n→∞ log |g(re iθ )| r ρ . (B6)
This describes the growth of g along a ray, {ξ : arg ξ = θ}. It is well known that h is continuous and that max

θ∈[0,2π] h(θ) = σ. (B7)
The first idea of GMR is that, with these notions, one can say that for large values of r, log g(re iθ ) ∼ h(θ)r ρ in some vague sense. This allows for quantizing the cancellation in function of the type σ and indicator function. According to [START_REF] Rubel | Entire and meromorphic functions[END_REF]Theorem 10.1], one can replace M g (r) by µ g (r) := max n |g n | r n in (B4), since log µ g (r) log M g (r) → 1, (r → ∞). Under certain additional assumptions [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF] and using (B6) one can estimate that d g (re iθ ) ∼ r ρ (σ g -h g (θ)). (B8) Equation (B8) is the main starting point in trying to find a way to evaluate the series (B1) without cancellation. One has to ensure that σ g -h g (θ) is as small as possible when evaluating along some ray or sector corresponding to θ. In the present case, θ = 0, since the evaluation is along the positive real axis.

The order, type and indicator of the probability function

This function, denoted by g, see Section III A, equation [START_REF] Alfano | Aerospace support to space situation awareness[END_REF], is an entire function of order ρ = 1, also called entire function of exponential type (EFET).

Proposition B.2. Let L g (λ) be the Laplace transform of g, given in equation [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF]. Let L g (λ) = Proof. This is proved using the correspondence between Borel and Laplace transforms, see for example [START_REF] Rubel | Entire and meromorphic functions[END_REF]Prop. 11.5,Cor. 11.5 and Prop. 11.7].

Finally, the indicator of g is obtained with Polya theorem, recalled below. (B10)

d. Reasons behind choosing to multiply by e pξ . The indicator function h g is showed in Figure 5 (a). One can see that in this case, d g (r) = σ g > 0 and thus, the sum is not optimally conditioned to be evaluated on the real axis. The GMR method suggests to multiply g by some preconditioner function ψ such that both ψ and f := ψg be very well conditioned for evaluation at the considered sector (positive real line in our case). If ψ is entire of order 1, we expect f (re iθ ) ∼ exp((h g (θ) + (B11)

Note that p -1 2σx 2 is always between the two other values, since by convention σ x ≥ σ y > 0.

Similarly, using Polya Theorem, Hence, from this analysis which follows closely the GMR method, the resulting series is optimally conditioned on the positive real axis for p ≥ Other choices are of course possible for ψ and currently, there is no established technique in the literature to assess which one is better. There are also "obvious bad choices": take for instance sin(ξ) which is entire of order 1 and indicator function h sin = ω |sin(θ)|, ω > 0. For real line evaluation, this will never give an optimal conditioning in zero.

h f (θ) = sup x∈ p-

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

2σ 2 y

 2 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

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.

Figure 3 :

 3 Figure 3: Number of terms in the series needed in Algorithm 2 for a requested accuracy δ ranging from 10 -1 to 10 -13 . Cases shown: Chan's Cases (left) and CSM cases (right) from Table1.

3 Chan 2 9 . 3 Chan 3 6 . 3 Chan 4 6 . 8 8 0 3 . 6 10 0 9 .

 3936368369 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 -3 Chan 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 Chan 7 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -8 6.443 × 10 -Chan 219 × 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 -Chan 656 × 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 -4 Chan 12 1.564 × 10 -9 1.564 × 10 -9 1.556 × 10 -9 1.564 × 10 -9 1.564 × 10 -9

CSM 3 7 .

 7 2004 × 10 -5 7.2000 × 10 -5 7.2000 × 10 -5 7.2003 × 10 -5 7.2003 × 10 -5

  are plotted in Figure4for ξ = 15. Observe that d g (ξ) 12 lost digits (taking log 10 ). Summing the first 101 terms with 10 digits precision in Maple, Digits:=10: add((-15.)^i/i!, i=0..100); gives -0.4847810247 • 10 -4 , while exp(-15) 3.06 • 10 -7 .

Figure 4 :

 4 Figure 4: Values of (-1) i 15 i i!

n→∞

  log log Mg(r) log r , where M g (r) := max |ξ|=r |g(ξ)|.

b.

  Relation between cancellation and type and indicator functions. Equation (B2) states that we are interested in minimizing the ratio d g (re iθ ) = log max n |g n | r n |g(re iθ )| .

  assuming σ x ≥ σ y > 0. Then, g(ξ) =

Theorem B. 3 .

 3 c. Borel Transform and Polya Theorem Let g be an EFET as above and the function L f (λ) = called the Inverse Borel Transform of g). This series converges outside the disk {|λ| ≤ σ g } and diverges inside this disk. The smallest convex compact set containing all its singularities is called the conjugate indicator diagram of g. We denote by k g (θ) the supporting function of this compact set. In general, the supporting function k(θ) of a set K ⊂ C is:k(θ) = sup z∈K {x cos θ + y sin θ} = sup z∈K { Re(ze -iθ )}, θ ∈ [0, 2π]. (B9)The following Theorem B.3 relates the growth of the EFET function g along a ray to the location of the singularities of its Inverse Borel transform. Polya[START_REF] Levin | Lectures on Entire Functions[END_REF] Chapter 9]. For every EFET g(ξ) the relation h g (θ) = k g (-θ) holds.Remark B.4. For our case, the singularities are S =

h

  ψ (θ))r). One possible choice is ψ(ξ) = e pξ for which the indicator function, h ψ = p cos(θ), θ ∈ [0, 2π] is given in Figure 5 (b) for p > 0. One can compute now the type and indicator of the preconditioned probability f , based on the Laplace transform of f given in Prop. III.1. The singularities are s0 = p, s1 = p -1 2σx 2 and s2 = p -1 2σy 2 , and hence σ f = max |p| , p -1 2σ y 2 .

σg 2 .

 2 We give for example the indicator function obtained when p = σg 2 in Figure5(c) and respectively p = σ g in Figure5 (d). For this problem, p = σ g is considered.

Figure 5 :

 5 Figure 5: Indicator functions of (a) probability g; (b) precondition ξ → e pξ ; (c) ξ → e σg 2 ξ g(ξ);

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[START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF]

	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, Algorithm 2) provide enough accuracy in most practical cases, as shown in Table 2. Except
	in Alfano's cases, they seem to be very effective due to their simplicity and may suffice in most
	practical cases.

Table 2 :

 2 Guaranteed correct digits obtained with formulas in Proposition III.6. The notation

e.g. 0.15

[START_REF] Chevillard | Multiple-Precision Evaluation of the Airy Ai Function with Reduced Cancellation[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.

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

	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

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.

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

and

.

Then:

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).

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

Proof of Proposition III.6

In general, for γ > 0,

and:

By multiplying the above expressions by e -pR 2 , one obtains:

) -1 p and:

Similarly, one gets: We consider a Taylor series of the form

where λ, α, κ > 0. For large ξ > 0, the computation in finite precision arithmetic of such sum is prone to cancellation. This is because the terms g i ξ i are first growing, before the series starts to converge when i κ ≥ αξ. When i κ αξ the terms g i ξ i usually get much larger than |g(ξ)|. So the leading bits cancel out while the lower-order bits which actually contribute to the first significant bits of the actual result get lost in the roundoff errors. So we are interested in minimizing the ratio )


Haiyang Ding

Mazen Alamir
email: mazen.alamir@grenoble-inp.fr

Francois Bonne

Ahmad Hably
email: ahmad.hably@grenoble-inp.fr

Patrick Bonnay

Distributed cooperative control framework of a cryogenic system

In this paper a recently proposed distributed control framework is applied to the model of a real-life cryogenic plant. The main advantage of the distributed control architecture is to allow for a modular design of the control algorithms. The paper shows that the cooperative nature of the solution enables an initial decentralized design to be improved by properly choosing the relative priority assignment between subsystems.

I. INTRODUCTION

Cryogenic systems are necessary to cool the supra-conducting devices that are used in many physical instruments including nuclear fusion reactors [START_REF] Henry | Analysis of the iter cryoplant operational modes[END_REF], [START_REF] Clavel | Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation[END_REF] and particle accelerators. Thanks to advances in system modeling in the last decade ( [START_REF] Bradu | Modeling, simulation and control of large scale cryogenic systems[END_REF], [START_REF] Bradu | A process and control simulator for large scale cryogenic plants[END_REF], [START_REF] Clavel | Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation[END_REF]), model-based control strategies have been used in order to avoid plants' oversize and to drastically reduce their energy consumption [START_REF] Bradu | Control optimization of a lhc 18kw cryoplant warm compression station using dynamic simulations[END_REF], [START_REF] Vinuela | Linear model-based predictive control of the lhc 1.8k cryogenic loop[END_REF], [START_REF] Clavel | Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation[END_REF], [START_REF] Bonne | Nonlinear observers of the thermal loads applied to the helium bath of a cryogenic joule-thompson cycle[END_REF]. While the above mentioned works are mainly based on centralized control schemes, the cryogenic system can be viewed as an interconnection of several coupled heterogeneous subsystems with local objectives. The cryogenic practitioners prefer rather modular design framework in which coupling-related issues and solutions are addressed by an additional layer that does not question the old existing widely assessed local controllers. This paper presents the performance of the application of the distributed control framework developed in [START_REF] Alamir | A novel distributed nmpc control structure for partially cooperative systems under limited information sharing[END_REF], [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] on a cryogenic system. This is realized by using the model of the experimental helium refrigerator facility located at CEA 1 -INAC 2 -SBT 3 , Grenoble, France. The paper is organized as follows: section II gives an overview of the cryogenic system. The model and the control objectives are presented in section III, while the cooperative distributed framework proposed in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] is briefly recalled in section IV. The simulation results and analysis are given in section V. Section VI concludes the paper.

II. OVERVIEW OF THE CRYOGENIC SYSTEM

The cryogenic system located at CEA-INAC-SBT, Grenoble is shown in Fig. 1. This plant offers a nominal capacity of 400 (watts W) at 4.5 (kelvin K) and serves as a testbed on which physical experiments are performed (testing cryogenic components, study of super-fluid helium, etc.) [START_REF] Clavel | Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation[END_REF]. Fig. 2 shows a block diagram of the cryogenic system consisting of a warm Compression system and a cold box (phase separator) in which the heating device emulates the devices to be cooled. Fig. 3 1 CEA: Commissariat à l'Energie Atomique et aux Energies Alternatives shows the ideal Claude thermodynamic cycle. As illustrated in Fig. 2, the Warm Compression System consists of a screw compressor with operational condition between 1.05 bar and 16 bar with a maximum flow rate of 72 g • s -1 and three control valves. The operational point of the refrigerator is set by the by-pass valve CV 956 . It is used to pass the excess flow rate towards the cold box. Control valves CV 952 and CV 953 serve respectively to supply or to evacuate the gas from the system via a helium gas drum. The heat load that disturbs the cryogenic system is represented by the resistance N CR 22 at the bottom of the diagram in Fig. 2. The cold box cools down the helium flow from 300 K to 4.5 K using the following equipment:

1) A liquid Nitrogen cooler followed by several counterflow heat exchangers. 2) A cold turbine expander, controlled by the CV 156 valve, which extracts work form the gas. 3) A Joule-Thomson expansion valve (CV 155 ). 4) A phase separator with the two-phase helium bath that is connected to the load.

Fig. 2: A block diagram of the cryogenic system consisting of a warm compression system and a cold box (phase separator) in which the heating device emulates the devices to be cooled

III. MODEL DESCRIPTION AND CONTROL OBJECTIVES

A. Model description

According to [START_REF] Majecki | Control-oriented modelling of the cold boc of the cryogenic station 400w@4.5k[END_REF], the system can be viewed as the interconnection of four subsystems (see Fig. 5): Joule-Thomson cycle (phase separator), Brayton cycle (Turbine), warm end, and warm zone. The cooperative distributed control is applied based on this subdivision of the cryogenic plant and these four subsystems are referred to by 1 , 2 , 3 and 4 in the sequel.

The subsystems are interconnected through their physical outputs listed in Table . I, where the subscripts T , P , M , H and C stand for temperature, pressure, mass of liquid flow, hot and cold respectively. For example, T H4 stands for the output hot temperature of subsystem 4 . Given the interconnection topology shown in Fig. 4, the outputs of subsystems are introduced. Note that y i→j indicates the output of subsystem i that affects subsystem j .

y 1→2 = [T C1 M H1 M C1 ] T , y 2→1 = [T H2 P H2 P C2 ] T y 2→3 = [T C2 M H2 M C2 ] T , y 3→2 = [T H3 P H3 P C3 ] T y 3→4 = [T C3 M H3 M C3 ] T , y 4→3 = [T H4 P H4 P C4 ] T
The structure of the differential equations governing each subsystem and the corresponding size of the state and control vectors are briefly sketched Below. For an exhaustive presentation including the matrix definitions, the reader can refer to [START_REF] Ding | On a partially cooperative distributed control framework with priority assignement[END_REF]. So one has the following subsystems:

Subsystem 1 -Phase separator ξ1 = A 1 ξ 1 + B u 1 B ν 1 u 1 ν + B 2→1 y 2→1 (1)
y 1→2 = C 1 ξ 1 + D 1   u 1 ν y 2→1   Fig. 3: Ideal Claude thermodynamic cycle
where ξ 1 ∈ R 17×1 is the state vector of subsystem 1 , u 1 ∈ R is the control input, i.e. the Joule-Thompson expansion valve and ν is the exogenous input (heat load).

Subsystem 2 -Brayton cycle/turbine ξ2 = A 2 ξ 2 + B u 2 u 2 + B 3→2 y 3→2 + B 1→2 y 1→2 (2)
y 2→1 y 2→3 = C 2→1 C 2→3 ξ 2 + D 2   u 2 y 3→2 y 1→2   with ξ 2 ∈ R 32×1 is the state vector. u 2 ∈ R is the control input of subsystem 2 i.e
. the turbine speed of the Brayton cycle in Fig. 5.

Subsystem 3 -Warm end ξ3 = A 3 ξ 3 + B u 3 u 3 + B 4→3 y 4→3 + B 2→3 y 2→3 (3)
y 3→2 y 3→4 = C 3→2 C 3→4 ξ 3 + D 3   u 3 y 4→3 y 2→3  
where ξ 3 ∈ R 50×1 is the state vector and u 3 ∈ R is the control input of the Nitrogen cooler shown in Fig. 5.

y 4→3 = C 4 ξ 4 + D 4 u 4 y 3→4
where ξ 4 ∈ R 2×1 is the state vector and u 4 ∈ R 2×1 is the control of the warm zone representing the control valves CV 952 and CV 953 in Fig. 2 which are used to supply or evacuate gas out of the system.

B. Control problem

The main control objective is to maintain the temperature of the liquid helium (T helium ) in the phase separator at 4.5 K despite disturbing heat load. Besides this primary task, several constraints regarding security and operational conditions have to be respected, namely:

• The output of the temperature of Brayton cycle 2 (T H2 in Fig. 4) must remain in the interval [7.8 K, level = 50%. The control architecture is defined in a way that each subsystem (with its local control loops) can operate when it is connected or disconnected to the remaining parts. This is particularly true for the compression zone. When the subsystems are interconnected, a cooperation has to take place where each controller's behavior is amended by signals coming only from its neighbors.

IV. THE DISTRIBUTED CONTROL ARCHITECTURE

In this section, the application of the cooperative distributed control framework [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] to the cryogenic plant is introduced. In order to do this, the dynamic equations ( 1), ( 2), ( 3) and ( 4) is transformed into the standard form used in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF].

A. Equation transformation

The standard form of the dynamic equations used in the cooperative control framework proposed in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] is given by:

ẋi = A i x i + B i u i + j∈I ← i (A j→i x j + B j→i u j ) (5)
A careful examination of equations ( 1), ( 2), ( 3) and ( 4) shows that they exhibit algebraic loop (for instance y 1→2 depends on y 2→1 to cite a single example). To eliminate these algebraic loops, a standard state extension technique is used in which additional auxiliary variables are defined that are associated to fast dynamics (1 sec is sufficiently fast for the problem at hand). The extended resulting model becomes:

ẋ1 = A 1 x 1 + B u 1 u 1 + B ν 1 ν + A 2→1 x 2 (6) ẋ2 = A 2 x 2 + B u 2 u 2 + A 1→2 x 1 + A 3→2 x 3 (7) ẋ3 = A 3 x 3 + B u 3 u 3 + A 2→3 x 2 + A 4→3 x 4 (8) ẋ4 = A 4 x 4 + B u 4 u 4 + A 3→4 x 3 (9)
in which u i , for i = 1, 2, 3, 4 and ν keep their meanings in (1)-(4), while x i , i = 1, 2, 3, 4 are the extended states with the following dimensions:

x 1 ∈ R 20×1 , x 2 ∈ R 38×1 , x 3 ∈ R 56×1 , x 4 ∈ R 5×1
The outputs y i can be expressed in terms of the extended states by y i = Ci x i (see [START_REF] Ding | On a partially cooperative distributed control framework with priority assignement[END_REF] for a complete description).

B. Control design 1) Local control design:

The local controller of each subsystem is designed at its nominal point as if it was alone, that is to say, using the following nominal equation:

ẋi = A i x i + B u i u i
The closed-loop feedback gain K i for subsystem i is calculated by LQR design minimizing the cost function:

J = ∞ 0 x T i (τ )Q i x i (τ ) + u T i (τ )R i u i (τ )dτ with Q i = p i • I + q i • C T i C i , R i = r i • I 4 .
The values of the parameters used in the nominal LQR design are given in Table.II 5 . Fig. 5: A schematic view of the cryogenic system in which the plant is viewed as the interconnection of 4 subsystems.

q i r i C i 1 1 10 3 1 s.t. T helium = C 1 x 1 2 0 1 10 -2 C 2 = Ĉ2 3 0 1 10 -2 C 3 = Ĉ3 4 0 1 10 2 C 4 = Ĉ4
The design is intentionally simple in order to show that the cooperative scheme is capable to recovery even with loosely locally designed loops. The resulting closed-loop feedback gain for nominal control is K i for subsystem i . Thus the control of subsystem i can be written as:

u i = -K i x i + v i ( 10
)
in which v i is the cooperative control term that is to be determined by the cooperative control algorithm.

2) Cooperative control design: The cooperative control proposed in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] involves the definition of priority coefficients π i j that describes how important subsystem j is viewed by subsystem i. Given the "chain-like" graph corresponding to the coupling expressed by ( 6)-( 9), the corresponding priority matrix takes the following form:

Π =           π 1 1 π 1 2 0 0 π 2 1 π 2 2 π 2 3 0 0 π 3 2 π 3 3 π 3 4 0 0 π 4 3 π 4 4           ; π j i ≥ 0 (11)
Recall that in the cooperative framework of [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF], when i is affected by j , the following quantity is sent by i to j

W i j := L j→i (x i ) = 2x i T P i (12)
in which P i is the corresponding solution to the Riccati equation when defining the LQR control of subsystem i . This information is crucial if j is willing to cooperate with i because W i j A j→i x j represents the additional term in the derivative of the Lyapunov function x T i P i x i of i that comes from the interaction with j . By doing so, j can incorporate all such terms coming from its neighbors in order to define a cooperative optimal control problem according to:

min vj v j 2 Rj + i∈I → j π i j W i j A j→i x j (13)
It is then shown in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] that by defining the extended state z j which gathers x j and all the terms {π i j W i j } i∈I → j coming from all the subsystems i that are affected by j (this is denoted by i ∈ I → j ), the cooperative optimal control problem viewed by j can be put in the following form (see [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] for the detailed expressions):

żi = A i z i + B i v i ( 14
)
while the cost function (13) rewritten in terms of the extended states z i in (14) becomes:

min vi v i T R i v j + z i T Q coop i z i (15)
for an appropriate definition of Q coop i see [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF]. By solving the quadratic problem (15), the discrete-time cooperative control term v i (k) = K v i z i (k) can be obtained and injected in [START_REF] Majecki | Control-oriented modelling of the cold boc of the cryogenic station 400w@4.5k[END_REF]. Note also that as far as linear networks of systems are concerned, [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] proposes a systematic procedure for priority assignment that guarantees the stability of the whole system under the cooperative control framework described above.

V. SIMULATIONS

Recall that the control objective is to reject heating disturbances that are applied in the bath in order to emulate the heat pulses coming from the supra-conducting devices. Fig. 6 shows the profile of the heating power injected to subsystem 1 in all the simulations illustrated in this section. A similar heating power profile has been used in [START_REF] Clavel | Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation[END_REF]. Besides the outputs of the subsystems, the evolution of nominal Lyapunov functions are also used to evaluate the performance of the cooperative distributed control. The nominal Lyapunov function of subsystem i is defined as:

V i = x i T P i x i (16)
in which P i is the solution to the corresponding Riccati equation when designing the nominal LQR regulator for subsystem i . Note that when the subsystem is at its operational point, the value of the nominal Lyapunov function is zero. Due to the fact that the four subsystems of the cryogenic plant are different, the comparison shown hereafter focuses on the so called normalized nominal Lyapunov functions (N-N Lyapunov function) that is defined by:

V s1 i = V s1 i V max i , V s2 i = V s2 i V max i
where s 1 and s 2 denote two different assignments of control framework parameters explained later, while V max i is the maximum value of the two nominal Lyapunov functions along their evolution for subsystem i in each set of comparison between the priority settings denoted by s 1 and s 2 , namely:

V max i = max k∈{1,••• ,Ns} max V s1 i (k) , max V s2 i (k)
where N s is the simulation length. In some simulations, the performances are compared by referring to the residual of the N-N Lyapunov functions defined by:

∆V i = V s1 i -V s2 i (17)
Note that when ∆V i is negative, then assignment s 1 is better than s 2 as far as the subsystem i is concerned.

The behavior of the cryogenic system under the decentralized control-loop (i.e. the cooperative control is deactivated by using v i = 0 for all i, denoted by "non-cooperative control") is compared with the distributed control with equal priority assignment Π equal (i.e π j i in (11) equals 1) and then with distributed cooperative control with high priority given to the warm compression subsystem Σ 4 (i.e. the priority matrix Π 4 is chosen with π j 4 = 2 (for j = 3, 4)). The results are shown in Fig. 7, the N-N Lyapunov function of cooperative control (V coop i ) is systematically less than that of non-cooperative control (V noncoop i ). The outputs of subsystems are illustrated in Figs. 8910. One can clearly notice that P C4 violates its constraints (P C4 ≤ 1.2 bar) (Fig. 8) in non-cooperative control while this constraint is fulfilled under cooperative control. By giving high priority to the warm zone ( 4 ) by Π 4 , the performances of all subsystems are improved when compared to the equal priority setting. So it can be concluded that the priority assignment Π 4 which emphasizes the importance of the warm zone allows improved effectiveness of the cooperative control framework.

VI. CONCLUSION

In this paper, the cooperative control proposed in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF] is applied to a cryogenic plant which can be viewed as the interconnection of four subsystems with individual tasks. Comparing to the non-cooperative control strategy in which only nominal controller is activated, the cooperative control with equal priority improves the performance of the cryogenic plant under non-modeled disturbance (heating power) at the phase separator. Moreover, it is found that by assigning high priority to the warm zone, the effectiveness of the cooperative control framework is further improved. Future work is currently in progress thanks to the French ANR-Cryogreen project in which nonlinear constrained local controllers are expected to be used. The objective of this extension is twofold: perform an explicit constraints handling and fully exploit the knowledge-based models that are sometimes nonlinear while avoiding the contamination of these nonlinearities on the overall system. This can be achieved by replacing the controlled nonlinear subsystems (under their nonlinear feedback) by appropriate resulting nonlinear model with the right response time and use the resulting linear models to assign the priority coefficient following the algorithms proposed in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF].

Acknowledgment. The authors are grateful to the financial support from the National French Research Agency (ANR) through the CRYOGREEN project but also through the former CHEOPS project.

APPENDIX

The details on equation transformation to eliminate the algebraic loops from the cryogenic plant are given in [START_REF] Ding | A distributed cooperative control scheme with optimal priority assignment and stability assessment[END_REF]. In this section only the definitions of the Ĉi are given:

Ĉ1 = [ C1→2 I 3×3 ] Ĉ2 = C2→1 I 3×3 O 3×3 C2→3 O 3×3 I 3×3 Ĉ3 = C3→2 I 3×3 O 3×3 C3→4 O 3×3 I 3×3 Ĉ4 = [ C4→3 I 3×3 ]
Note that the notation O a×b and I a×b represent respectively the zero matrix and the identity matrix of size (a × b). Fig. 7: Evolution of the N-N Lyapunov functions of the subsystems under non-cooperative control (solid red), cooperative control with equal priority (dashed blue), and cooperative control with priority for warm zone. In general, the N-N Lyapunov function of cooperative control is less than that of the non-cooperative control. Fig. 8: P C4 and P H4 under decentralized control (solid red) cooperative control with equal priority (dashed blue) and cooperative control with priority for the warm zone (violet). In the non-cooperative control, P C4 goes beyond its upper limit (1.2 bar). In addition, the performance is improved with cooperative control. When the warm zone is given a high priority P C4 and P H4 are better regulated by cooperative control.

Fig. 1 :

 1 Fig. 1: Photos of the cryogenic plant of CEA-INAC-SBT, Grenoble. (a) The compressor of the warm compression system. (b) Global view of the cold box. (c) Details of the cold box.

Fig. 4 :

 4 Fig.4: The interconnection terms between the subsystems of the cryogenic plant.

Fig. 6 :

 6 Fig.6: Heating power profile injected into the cryogenic system.

Fig. 9 :

 9 Fig.9: P C2 and P H2 achieve improved performances although high priority is assigned to the warm zone.

Fig. 10 :

 10 Fig.10: The outputs of subsystem 3 with high priority to the warm zone.

TABLE I :

 I Outputs of the subsystems of the cryogenic plantBrayton cycleT H 2 , P H 2 , P C 2 , T C 2 , M H 2 , M C 2 , P H 3 , P C 3 , T C 3 , M H 3 , M C 3

	Symbol Subsystem name	Output
	3	warm end	T H 3

1 phase separator T C 1 , M H 1 , M C 1 2 4

warm zone T H 4 , P H 4 , P C 4

TABLE II :

 II Values of the parameters used in the nominal LQR design

	Subsystem p i

I stands for the identity matrix with a suitable dimension

See the appendix for the definitions of Ĉi