D Arzelier
email: arzelier@laas.fr

C Louembet
email: clouembet@laas.fr

A Rondepierre
email: aude.rondepierre@insa-toulouse.fr

M Kara-Zaitri

M Kara-Zaïtri

A New Mixed Iterative Algorithm to Solve the Fuel-Optimal Linear Impulsive Rendezvous Problem

Keywords: Orbital rendezvous, fuel optimal space trajectories, primer vector theory, impulsive maneuvers, linear equations of motion

The optimal fuel impulsive time-fixed rendezvous problem is reviewed. In a linear setting, it may be reformulated as a non convex polynomial optimization problem for a pre-specified fixed number of velocity increments. Relying on variational results previously published in the literature, an improved mixed iterative algorithm is defined to address the issue of optimization over the number of impulses. Revisiting the primer vector theory, it combines variational tests with sophisticated numerical tools from algebraic geometry to solve polynomial necessary and sufficient conditions of optimality. Numerical examples under circular and elliptic assumptions show that this algorithm is efficient and can be integrated into a rendezvous planning tool.

Introduction

Given the increasing need for satellite servicing in current and future space programs developed in conjunction with rendezvous missions for the International Space Station (ISS), the interest of most space agencies in developing adequate rendezvous mission planning tools has been rapidly rising. In particular, new challenges have appeared relating to the synthesis of guidance schemes capable of achieving autonomous far range rendezvous on highly elliptical orbits while preserving optimality in terms of fuel consumption. Strictly speaking, the space far range rendezvous maneuver is an orbital transfer between a passive target and an actuated spacecraft called the chaser, within a fixed or floating time period. In this paper, we mainly focus on the so-called time-fixed fuel optimal rendezvous problem in a linearized gravitational field for which a renewed interest has been witnessed in the literature [START_REF] Hughes | A comparison of trajectory optimization methods for the impulsive minimum fuel rendezvous problem[END_REF], [START_REF] Larsson | Fuel Efficient Relative Orbit Control Strategies for Formation Flying and Rendezvous within PRISMA[END_REF], [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF], [START_REF] Delpech | Preliminary results of the vision based rendezvous and formation flying experiments performed during the PRISMA extended mission[END_REF]. Because of the constraints of on-board guidance algorithms, numerical solutions based on linear relative motion are particularly appealing. With respect to the numerical solution, direct methods based on discretizing the original problem and converting it into a linear programming problem may be used as in [START_REF] Larsson | Fuel Efficient Relative Orbit Control Strategies for Formation Flying and Rendezvous within PRISMA[END_REF]. Indirect approaches based on the solution of optimality conditions derived from Pontryagin's Maximum Principle, leading to the development of the so-called primer vector theory presented in [START_REF] Lawden | Optimal trajectories for space navigation[END_REF], have also been an avenue of research in numerous studies [START_REF] Prussing | Optimal multiple-impulse orbital rendezvous[END_REF], [START_REF] Prussing | Illustration of the primer vector in time-fixed orbit transfer[END_REF], [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], [START_REF] Carter | Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion[END_REF]. As they only focus on fixed number and location of impulses, these approaches fail to optimize trajectory planning in terms of the number of impulsive maneuvers. Neither can they ensure the global optimality of the result provided for a fixed number of impulses.

To optimize the number of impulses as well as their specific application times, an iterative algorithm based on the calculus of variations originally developed by Lion and Handelsman [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF] has been designed [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF][START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF]. The main drawback however is due to the possible non smoothness and sub-optimality of the resultant trajectory of the primer vector norm due to the local nature of the proposed approach. To overcome this difficulty, a Davidon-Fletcher-Powell penalty minimization step is proposed in order to move the impulses and achieve a smooth trajectory as detailed in [START_REF] Prussing | Optimal multiple-impulse time-fixed rendezvous between circular orbits[END_REF] or [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF]. The paper's contribution is to revisit the iterative algorithm of Lion-Handelsman based on the calculus of variations and to take advantage of the polynomial nature of the underlying necessary conditions to circumvent the necessity to resort to local optimization schemes. In particular, this new iterative algorithm combines the algebraic formulation of Carter's necessary and sufficient conditions, supporting the use of powerful numerical tools from the algebraic geometry field, and improving variational tests as derived in [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF]. This results in a mixed different iterative strategy bypassing the local optimal search step and cusp occurrence. Circular and elliptic Keplerian reference orbits may be considered indiferrently.

In the first section of this paper, the framework of the minimum-fuel fixed-time rendezvous problem is presented and necessary and sufficient conditions of optimality are recalled. Relative dynamics motion for rendezvous are the well-known Tschauner-Hempel equations [START_REF] Tschauner | Rendevous zu einem in elliptischer Bahn umlaufenden Ziel[END_REF] and the transition matrix of Yamanaka-Ankersen [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF]. The results of [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF] are recalled and the mixed iterative algorithm is presented. For comparison's sake, the efficiency of the proposed algorithms is illustrated with four different numerical examples. Three academic examples taken from Carter's reference [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] are first studied in details. Also reviewed is one realistic scenario based on PRISMA which is a "technology in-orbit testbed mission" demonstrating formation flight [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF].

The Time-Fixed Optimal Rendezvous Problem

Linear impulsive time-fixed optimal rendezvous problem

This paper focuses on the fixed-time minimum-fuel rendezvous between close orbits of an active (actuated) spacecraft called the chaser with a passive target spacecraft assuming a linear impulsive setting and a Keplerian relative motion as it is defined in the references [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], [START_REF] Carter | Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion[END_REF]. The impulsive approximation for the thrust means that instantaneous velocity increments are applied to the chaser whereas its position is continuous. If the relative equations of motion of the chaser are supposed to be linear and under the previous Keplerian assumptions, it is shown in the references [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], [START_REF] Carter | Linearized impulsive rendezvous problem[END_REF], [START_REF] Carter | Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion[END_REF] that the considered minimum-fuel rendezvous problem may be reformulated as the following optimization problem: min

N,ν i ,∆ v i ,β (ν i ) J = N ∑ i=1 ∆ v i s.t. u f = N ∑ i=1 φ -1 (ν i )B(ν i )∆ v i β (ν i ) = N ∑ i=1 R(ν i )∆ v i β (ν i ) β (ν i ) = 1 ∆ v i ≥ 0 (1)
where φ (ν) is the fundamental matrix associated to the linearized relative free motion and Φ(ν, ν 1 ) = φ (ν)φ -1 (ν 1 ) denotes therefore the transition matrix of the linearized relative free motion. Note that the true anomaly ν has been chosen as the independent variable throughout in the paper. ν 0 and ν f denote respectively the initial and final values of the true anomaly during the rendezvous.

u f = φ -1 (ν f )X f -φ -1 (ν 1 )X 0 = 0 where the state vector X f = r T f v T f T
at ν f and the state vector

X 0 = r T 0 v T 0 T
at ν 0 are composed of the relative positions and relative velocities vectors. The optimization decision variables are the number of impulses N, the sequence of thrust times

{ν i } i=1,••• ,N , the sequence of thrust magnitudes {∆ v i } i=1,••• ,N and of thrust directions {β (ν i )} i=1,••• ,N .
Due to the lack of a priori information about the optimal number of impulses to be considered, problem [START_REF] Hughes | A comparison of trajectory optimization methods for the impulsive minimum fuel rendezvous problem[END_REF] is very hard to solve from both theoretical and numerical points of view. Therefore, the associated fixed-time minimum-fuel rendezvous problem for a fixed number N of impulses has been considered in the literature mainly via geometric methods near circular [START_REF] Prussing | Optimal multiple-impulse orbital rendezvous[END_REF], [START_REF] Prussing | Optimal multiple-impulse time-fixed rendezvous between circular orbits[END_REF], [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF] or elliptic [START_REF] Carter | Linearized impulsive rendezvous problem[END_REF] orbits. These results are mainly based on the derivation of optimality conditions for the problem (1) when N is fixed a priori.

Carter's necessary and sufficient conditions for a fixed number of impulses

When the number of impulses is not a part of the optimization process and is fixed a priori to N, problem (1) may be considered as the joint optimal selection of N velocity increments ∆V (ν i ) = ∆ v i β (ν i ) and N times ν i of maneuvers. Applying a Lagrange multiplier rule for the problem (1) as in [START_REF] Carter | Linearized impulsive rendezvous problem[END_REF], one can derive necessary conditions of optimality (2) to [START_REF] Prussing | Optimal multiple-impulse orbital rendezvous[END_REF] in terms of the Lagrange multiplier vector λ ∈ R n as is recalled in theorem 2.1 below. Prussing has first shown in [START_REF] Prussing | Optimal impulsive linear systems: Sufficient conditions and maximum number of impulses[END_REF] that these conditions are also sufficient in the case of linear relative motion with the strengthening semi-infinite constraint (9) that should be fulfilled on the continuum [ν 0 , ν f ] and is expressed in terms of the transition matrix R(ν).

Theorem 2.1 [9]

(ν 1 , ..., ν N , ∆ v 1 , ..., ∆ v N , β (ν 1 ), ..., β (ν N )) is the optimal solution of problem (1) if and only if there exists a non-zero vector λ ∈ R m , m = dim(φ ) that verifies the necessary and sufficient conditions:

∆ v i = 0 or β (ν i ) = R T (ν i )λ , ∀ i = 1, • • • , N (2)
∆ v i = 0 or λ T R(ν i )R(ν i ) T λ = 1, ∀ i = 1, • • • , N (3)
∆ v i = 0 or ν 1 = ν 0 or ν N = ν f or λ T dR(ν i ) dν R(ν i ) T λ = 0, ∀ i = 1, • • • , N (4) N ∑ i=1 R(ν i )R T (ν i ) λ ∆ v i = u f (5) ∆ v i ≥ 0, ∀ i = 1, • • • , N (6)
N ∑ i=1 ∆ v i = u T f λ > 0 (7)
u T f λ is the minimum of the set defined as : {λ ∈ R m : (2) -(7) are verified} (8)

λ v (ν) ≤ 1 ∀ ν ∈ [ν 1 , ν N ] (9)
where λ v (ν) = R(ν) T λ denotes the so-called primer vector. Note that conditions [START_REF] Prussing | Illustration of the primer vector in time-fixed orbit transfer[END_REF] and ( 8) may be easily derived from the previous ones. These results derive directly from the seminal work of [START_REF] Lawden | Optimal trajectories for space navigation[END_REF] in the early sixties and form an alternative formulation to the primer vector theory. The primer vector is nothing but the co-state vector associated with the relative velocity and is related to the Lagrange multiplier vector through the preceding relation.

A numerical solution of optimality conditions (2) to [START_REF] Carter | Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion[END_REF] in the unknowns

λ ∈ R m , {ν i } i=1,••• ,N , {β (ν i )} i=1,••• ,N , {∆ v i } i=1,••• ,N
is still hard to find for a fixed number of impulses N due to the non convex and transcendental nature of these polynomial equalities and inequalities. However, it is shown in [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] that the problem may be tackled by using an adequate dynamic gridding strategy and polynomial optimization. Still, the optimal number of impulses N * for a particular rendezvous problem is in general unknown and only a bound N * ≤ N max is available [START_REF] Neustadt | A general theory of minimum-fuel space trajectories[END_REF]. N max = 2 for out-of-plane rendezvous, N max = 4 for in-plane rendezvous while N max = 6 for a general three-dimensional rendezvous problem. Apart the gridding step, the necessity to try every case for 2 ≤ N ≤ N max in [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] appears be very time-consuming and a more direct approach to solve problem ( 1) is now proposed.

3 Optimizing over the Number of Impulses

Using Lion & Handelsman results on multi-impulse trajectories

In [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF], a method is proposed to take advantage of the primer vector theory developed by Lawden in order to improve non optimal trajectories by adding or shifting impulses. The calculus of variations is used to find conditions on the norm of the primer vector for an additional impulse and on the derivative of this norm for initial and/or final coastings. The method is mainly based on derivation of the so-called variational adjoint equation resulting from the variation of the cost function. Later, Jezewski [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF], [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF] developed a numerical algorithm combining Lion-Handelsman's conditions with a modified gradient search approach in a linear model setting. The additional local optimization procedure is used to find the optimal position and modulus of the additional impulse so as to avoid a resulting cusp for the norm of the primer vector as reported in [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF], [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF]. In this section, the conditions of Lion and Handelsman are recalled and a different Heuristic iterative procedure avoiding local optimal search step and cusp occurrence is proposed. It is worth noticing that the extension of these conditions for elliptical reference orbit (Tschauner-Hempel [START_REF] Tschauner | Rendevous zu einem in elliptischer Bahn umlaufenden Ziel[END_REF] dynamical relative model and Yamanaka-Ankersen transition matrix [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF]) from [START_REF] Kara-Zaitri | Mixed iterative algorithm for solving impulsive time-fixed rendezvous problem[END_REF] is used here.

Additional interior impulse condition

Perturbing a reference initial two-impulse trajectory and adding an interior impulse at ν m , the differential cost can be expressed as:

δ J = ∆ v m (1 -λ v (ν m ) T β (ν m )) = ∆ v m (1 -λ T R(ν m )β (ν m )) (10)
From [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF], it is easy to conclude that δ J < 0 when λ v (ν m ) > 1 and that a maximum decrease in cost is obtained when:

ν m = arg max ν∈[ν 0 ,ν f ] λ v (ν) = arg max ν∈[ν 0 ,ν f ] R T (ν)λ (11)

Additional coasting period conditions

For an additional initial coasting period of duration dν 1 , the cost variation is given by:

δ J = -∆ v 1 dλ v dν T (ν 1 )λ v (ν 1 )dν 1 = -∆ v 1 λ T dR dν (ν 1 )R T (ν 1 )λ dν 1 = ∆ v 2 1 d λ v (ν) dν |ν=ν 1 dν 1 ( 12
)
This condition means that adding an initial coasting arc of dν 1 > 0 duration may improve the cost if λv (ν 1 ) T λ v (ν 1 ) > 0, i.e., the right derivative of the primer vector norm at ν 1 is positive. Similarly, for a final coasting arc of duration dν f , we get:

δ J = -∆ v f dλ v dν (ν f ) T λ v (ν f )dν f = -∆ v f λ T dR dν (ν f )R T (ν f )λ dν f = -∆ v 2 f d λ v (ν) dν |ν=ν f dν f (13)
A final coast of dν f < 0 duration will improve the cost when the left derivative of the primer vector norm at ν f is negative, i.e., when

dλ v dν (ν f ) T λ v (ν f ).

A mixed iterative algorithm

These conditions may be used jointly to reduce the cost of a reference non optimal two-impulse trajectory but can also be generalized to multi-impulse trajectories. Consider the four-impulse trajectory of Fig. 1a. Overall cost can be reduced by considering the coasting arc [ν i , ν i+1 ] in different ways:

-Adding a new impulse at ν m , -Adding an initial coasting arc by shifting

ν i toward ν m , if d λ dν (ν i ) > 0,
-Adding a final coasting arc by shifting

ν i+1 toward ν m , if d λ dν (ν i+1 ) < 0, -Replacing ν i and ν i+1 by ν m , if d λ dν (ν i ) > 0 and d λ dν (ν i+1 ) < 0.
This is equivalent to adding an initial and a final coasting arc on [ν i , ν i+1 ] and an impulse at ν m . As noted in [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF] and in [START_REF] Guzmán | Primer vector optimization: Survey of theory, new analysis and applications[END_REF], computation of the mid-impulse might nevertheless result in a non optimal trajectory not verifying the optimality conditions of Lawden and condition [START_REF] Carter | Necessary and sufficient conditions for optimal impulsive rendezvous with linear equations of motion[END_REF], particularly in the case of occurrence of a cusp at ν m as illustrated in Fig. 1b. A particular strategy combining Lion-Handelsman's conditions and local direct optimization based on the Davidon-Fletcher-Powell penalty method in [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF] or based on BFGS method in [START_REF] Guzmán | Primer vector optimization: Survey of theory, new analysis and applications[END_REF] has been proposed to optimize the resulting three-impulse trajectory. The objective of this section is to propose an alternative to this complicated procedure by developing a mixed iterative algorithm taking advantage of the algebraic formulation of Carter's optimality conditions and of the Lion-Handelsman's conditions. Starting from a non optimal two-impulse trajectory, successive admissible improved trajectories will be iteratively built by:

|λ v | ν ν 1 ν i ν i+1 ν f ν m 1 (a)
-Adding impulse at ν m if the impulse number does not exceed the upper bound N max ; -Moving the proximal impulse to ν m if an impulse cannot be added due to N = N max ; -Merging two impulses at ν m if there is no proximal impulse and if an impulse cannot be added due to N = N max .

where ν m is defined by [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF]. The logic of the proposed heuristic algorithm is depicted at Figure 2 where d

p i = d λ dν (ν i ) and d p i+1 = d λ dν (ν i+1
). The set T imp denotes the current set impulses and the new impulse ν m is always added to it.

d p i > 0 d p i+1 < 0 ν i = ν 0 ν i+1 = ν f νm = νi + νi+1 2 ν i+1 = ν m and T imp = T imp -{ν i+1 } ν i+1 = ν m and T imp = T imp -{ν i+1 } ν i = ν m and T imp = T imp -{ν i } ν i = ν m and T imp = T imp -{ν i } ν i = ν i+1 = ν m and T imp = T imp -{ν i , ν i+1 } ν i = ν i+1 = ν m and T imp = T imp -{ν i , ν i+1 } ν i = ν m or ν i+1 = ν i and T imp = T imp -{ν i or ν i+1 }
Fig. 2: Heuristic for the iterative mixed algorithm This heuristic procedure relies on basic principles that are used to make the successive sequences of maneuvers monotonically converging to an optimal solution.

-The first and final maneuvers time defined respectively at ν 0 and ν f cannot be moved in the process. Therefore, optimal solutions consisting of an initial and/or final coasting period cannot be found by this algorithm.

-A new impulse is always added at ν m . In case 8 for an actual number of impulses N < N max , it increases the number of impulses while when N = N max it reduces the number of impulses in cases 3 and 7 and it does not change this number for cases 1, 2, 4, 5 and 6.

-Every move of an actual maneuver time to ν m is chosen to be the proximal of ν m (in particular for case 4) except in case 7 which is specific since moves of impulses ν i and ν i+1 are not required to improve the trajectory. In this case, the idea is rather to re-initialize the iterative process with a new three-impulse trajectory. In practice, this case has never occurred on all the tested numerical examples.

The systematic convergence of the algorithm for any rendezvous is not analytically proved but no such case has been reported in the different numerical tests performed so far.

The algorithm may now be described in details. It is mainly composed of two stages: one initialization step solving a twoimpulse rendezvous problem and the iterative procedure building the final plan of maneuvers. ε cond and ε λ are respectively precision values on the conditioning number of the transition sub-matrix Φ 12 and on the maximum of the norm of the primer vector. Typically, ε cond = 10 6 while ε λ = 10 -6 .

Initialization step:

a. Solve the two-impulse problem:

1-Initialize

T imp = {ν 1 , ν 2 } = ν 0 , ν f (14)
2-Compute the transition matrices

Φ(ν f , ν 1 ) =     Φ 11 (ν f , ν 1 ) Φ 12 (ν f , ν 1 ) Φ 21 (ν f , ν 1 ) Φ 22 (ν f , ν 1 )     , Φ # (ν f , ν 1 ) = Φ -T (ν f , ν 1 ) =     Φ # 11 (ν f , ν 1 ) Φ # 12 (ν f , ν 1 ) Φ # 21 (ν f , ν 1 ) Φ # 22 (ν f , ν 1 )     (15)
If cond Φ 12 (ν f , ν 1 ) < ε cond Then ∆V (ν 1 ) = Φ -1 12 (ν f , ν 1 ) r f -Φ 11 (ν f , ν 1 )r 1 -Φ 12 (ν f , ν 1 )v 1 ∆V (ν f ) = v f + Φ 22 (ν f , ν 1 )Φ -1 12 (ν f , ν 1 )Φ 11 (ν f , ν 1 ) -Φ 21 (ν f , ν 1 ) r 1 -Φ 22 (ν f , ν 1 )Φ -1 12 (ν f , ν 1 )r f (16) Else Solve polynomial system w.r.t. (λ , {∆ v i } i=1, f ) λ T R(ν i )R(ν i ) T λ = 1, ∀ ν i ∈ T imp ∑ ν i ∈T imp R(ν i )R T (ν i ) λ ∆ v i = u f , ∀ ν i ∈ T imp , i = 1, f ∆ v i ≥ 0, i = 1, f (17)
Choose the minimum-fuel solution:

λ init = arg min λ u T f λ (18)
Compute impulses:

β (ν i ) = R(ν i ) T λ init , ∀ ν i ∈ T imp , i = 1, f ∆V (ν i ) = ∆ v i β (ν i ), ∀ ν i ∈ T imp , i = 1, f (19) b. Propagate primer vector λ v (ν) on a grid Π = {ν 1 , • • • ν f }: If cond Φ # 21 (ν f , ν 1 ) < ε cond Then λ v (ν 1 ) = ∆V (ν 1 ) ∆ v 1 λ v (ν f ) = ∆V (ν f ) ∆ v f λ v (ν) = Φ # 21 (ν, ν 1 )Φ #-1 21 (ν f , ν 1 ) λ v (ν f ) -Φ # 22 (ν f , ν 1 )λ v (ν 1 ) + Φ # 22 (ν, ν 1 )λ v (ν 1 ) (20)

Else

Compute R(ν) on grid Π and Propagate λ v (ν) via

λ v (ν) = R T (ν)λ (21)
c. Compute

λ vm = max ν∈[ν 0 ,ν f ] λ v (ν) and ν m = arg max ν∈[ν 0 ,ν f ] λ v (ν) d. If λ vm -1 < ε λ Then stop.
The two-impulse trajectory is optimal.

Else start iterative procedure.

ε 1λ and ε 2λ are two different precision values used in the iterative procedure to test the maximum value of the norm of the primer vector with respect to 1. ε cost is a parameter used to check the evolution of the cost during the iterative process.

Iterative procedure:

While (λ vm -1 ≥ ε 1λ ) and diff cost = |u T f λ iter -u T f λ iter-1 | > ε cost or λ vm -1 ≥ ε 2λ ) a. iter← iter+1 ; Choose ν a , ν b ∈ T imp such that (ν a < ν m < ν b ) and (ν a = ν i , ν b = ν i+1 ) b. Modify T imp = T imp ∪ {ν m }. If dim(T imp ) > N max Then Modify T imp as (1) If dλ v (ν a ) dν T λ v (ν a ) > 0 and dλ v (ν b ) dν T λ v (ν b ) < 0 Then (i) If ν a = ν 1 Then ν b = ν m and T imp = T imp -{ν b } (22) (ii) Else If ν b = ν f Then ν a = ν m and T imp = T imp -{ν a } (23)
(iii) Else

ν a = ν b = ν m and T imp = T imp -{ν a , ν b } if ν m = (ν a + ν b )/2 (24) ν a or ν b = ν m and T imp = T imp -{ν a or ν b } if |ν m -ν a | < |ν m -ν b | or |ν m -ν b | < |ν m -ν a | (25)
(2)

If dλ v (ν a ) dν T λ v (ν a ) < 0 or dλ v (ν b ) dν T λ v (ν b ) > 0 (i) If dλ v (ν a ) dν T λ v (ν a ) > 0 Then ν a = ν m and T imp = T imp -{ν a } (26) (ii) If dλ v (ν b ) dν T λ v (ν b ) < 0 Then ν b = ν m and T imp = T imp -{ν b } (27) (iii) If dλ v (ν b ) dν T λ v (ν b ) > 0 and dλ v (ν a ) dν T λ v (ν a ) < 0 ν a = ν b = ν m and T imp = T imp -{ν a , ν b } (28) c. Solve polynomial system w.r.t. (λ , {∆ v i } i=1, f ) λ T R(ν i )R(ν i ) T λ = 1, ∀ ν i ∈ T imp ∑ ν i ∈T imp R(ν i )R T (ν i ) λ ∆ v i = u f , ∀ ν i ∈ T imp , i = 1, • • • , N ∆ v i ≥ 0, i = 1, • • • , N (29)
d. Choose the minimum-fuel solution:

λ iter = arg min λ u T f λ (30)
e. Compute impulses:

β (ν i ) = R(ν i ) T λ iter , ∀ ν i ∈ T imp , i = 1, • • • , N ∆V (ν i ) = ∆ v i β (ν i ), ∀ ν i ∈ T imp , i = 1, • • • , N (31) f. Compute cost difference diff cost = |u T f λ iter -u T f λ iter-1 | (32) g. Compute R(ν) on grid Π and Propagate λ v (ν) via λ v (ν) = R T (ν)λ iter (33) h. Compute λ vm = max ν∈[ν 1 ,ν f ] λ v (ν) and ν m = arg max ν∈[ν 1 ,ν f ] λ v (ν)
Repeat Iterative Procedure until λ vm ≤ 1

Initialization stage and step c of iterative procedure require solving a system of polynomial equations ((17) in the first case if the transition matrix Φ(ν f , ν 1 ) is ill-conditioned and (29) in the iterative procedure) with respect to λ and ∆ v i . The set of solutions to the two first set of equations of this polynomial system is composed of 8 couples of solutions (λ

, {∆ v i } i=1,••• , f ).
Among this set of solutions, only those corresponding to a positive magnitude are kept to compute the minimum-fuel solution.

Note that regular algebraic tools for finding all real solutions of multivariate polynomial equations based on formal Gröbner basis computation may fail due to highly complex equations. Here, homotopy continuation methods have been used [START_REF] Verschelde | Homotopy methods for solving polynomial systems[END_REF]. In particular, the free software package PHCPack developed by Jan Verschelde [START_REF] Verschelde | Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation[END_REF] is used to solve the system of polynomial equations at each iteration at step c of the iterative procedure and at the initialization step if necessary. The efficiency of this algorithm is now demonstrated on several different examples.

Applications and Numerical Examples

In this section, numerical results obtained from the mixed iterative algorithm are compared with previous ones published in the literature [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] and on a more realistic case borrowed from the PRISMA test bed [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF]. When the mixed iterative algorithm converges to an N it solution, then the PRDV algorithm from [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] is used to certify optimality of this solution for this fixed number of impulses. Only coplanar elliptic rendezvous problems based on the Yamanaka-Ankersen transition matrix [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF] are considered for numerical illustration of the results proposed. Note also that the algorithm has been successfully applied to the highly elliptic rendezvous mission SIMBOL-X in [START_REF] Kara-Zaitri | Mixed iterative algorithm for solving impulsive time-fixed rendezvous problem[END_REF]. In this case and under Keplerian assumptions, the bound of Neustadt [START_REF] Neustadt | A general theory of minimum-fuel space trajectories[END_REF] on the optimal number of impulses is 4 and therefore N max = 4 in the following. Note that for the first three particular cases in which eccentricity e = 0, the Yamanaka-Ankersen transition matrix reduces to the Hill-Clohessy-Wiltshire transition matrix [START_REF] Clohessy | Terminal guidance system for satellite rendezvous[END_REF]. Finally, all numerical examples are processed using PHCpack 2.3.52 [START_REF] Verschelde | Homotopy methods for solving polynomial systems[END_REF], [START_REF] Verschelde | Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation[END_REF] under Matlab 2010b c running on an Intel Core(TM) i7 X920 2.0GHz system with 8GB ram.

Case study 1

The first numerical example is defined by the third case first considered in the reference [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] and later in [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] via polynomial optimization. Data given in Table 1 define a circle-to-circle rendezvous problem where the chaser is one unit before the target with the same velocity and has to rendezvous with the target with a positive vertical velocity after one orbital period.

Eccentricity e = 0 ν 0 0 rad X T 0 [ 1 0 0 0 ] ν f 2π rad X T f [ 0 0 0 0.427 ]

Table 1: Data for Carter's third example

The mixed iterative algorithm reaches the optimal four-impulse solution within 16 seconds after 6 iterations with an initial and a final impulse. For comparison's sake, the Table 2 recalls the results given in the references [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] and [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] alongside the ones given by the mixed iterative algorithm and obtained from the application of the procedure given in [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF]. Remember that Carter's results are obtained using a priori pre-assigned thrust times. The results of the mixed iterative algorithm may be certified via the comparison with those of [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] and by applying the quadratic-based method proposed in [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF] to build 4-impulse optimal solutions for circular rendezvous.

Carter [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] Mixed iterative algorithm PRDV algorithm [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] Analytic procedure [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF] Note that the results of PRDV algorithm are rather optimistic in terms of consumption due to numerical approximations in the computation process leading to a solution that does not respect exactly the symmetry requirements of four-impulse optimal solutions [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF]. Norm of the primer vector, thrust directions, primer vector locus, cost evolution during the iterative process are depicted on Figure 3.

λ v (ν) β(ν)
Fig. 3: λ v (ν) , primer vector locus, L 1 /l 2 cost for the third case of Carter's example [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] Figures 4 to 5 show the iterative construction of the optimal solution following the proposed strategy. Note that the norm of the primer vector λ v (ν) is a smooth function of the anomaly at every step of the process, thereby overcoming the main drawback of the usual iterative procedure originally proposed by Lion-Handelsman [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF] and developed lately in [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF] and [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF].

Red impulses give the localization of the thrusts effectively computed at the preceding step while the green impulse gives the localization of the resulting impulse from either addition of a new impulse or move of a red impulse. Fig. 4: Details of iterations 1-3 for the third case of Carter [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] The whole process is formed of two separate stages as by is easily deduced from the analysis of these details. The first stage consists mainly in a successive addition of new impulses from the initial two-impulse solution to a non-optimal four-impulse solution. Once the upper-bound of Neustadt is reached, the algorithm moves the interior impulses to find their optimal location. The second numerical illustration is still a coplanar circle-to-circle rendezvous that should be completed within one orbital period and that is borrowed from the reference [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]. The boundary conditions of the rendezvous are modified so that the chaser is one unit above the target with the same initial and final velocities. All characteristics are summarized in Table 3.

Eccentricity e = 0 ν 0 0 rad X T 0 = [r T 0 v T 0 ] [ 0 1 0 0 ] ν f 2π rad X T f = [r T f v ′ f ] [ 0 0 0 0 ] N max 4
Table 3: Carter's second example characteristics

The resulting optimal trajectory computed with the mixed iterative algorithm is a three-impulse rendezvous as conjectured by Carter. The primer vector trajectory plot depicted in Fig. 7 confirms the optimality of the solution for a fixed number of impulses equal to N = 3.

-1.5 1.5

λ vx λ vz Unit circle λ v (t) β (t i ) 0 0 0.5 0.5 1.0 -0.5
-0.5 -1.0 Fig. 7: Primer vector in-plane trajectory for Carter's case study 2 [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] The optimization process requires 9 iterations during 25 s. These ten steps are detailed in Figures 8 through 10, where each iteration is associated with one particular T imp update case. As previously, the whole process may be divided into two different stages. The first one is similar to the first example while the second one is made off two alternate operations: Addition of a fourth interior impulse to the previous three-impulse non optimal solution and merging of two interior impulses to get a better three impulse solution. It is important to note that there is no admissible extreme two-impulse solution since the following linear system corresponding to the two-point boundary value problem has no solution. The results presented in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] are clearly not optimal with respect to the choice of the date of application of the interior velocity increment. In [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], this location has been chosen a priori and it obviously results in a non optimal consumption as summarized in Table 4 where the results of the mixed iterative algorithm are compared to those presented by Carter in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] with 12% fuel consumption improvement. Additionally, Table 4 lists the optimality certificate furnished by the PRDV algorithm [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF].

            0 2 6π 2 0 1 0 1 2 0 2 0 -1 0 -1 0                         ∆V x (0) ∆V z (0) ∆V x (2π) ∆V z (2π)             =             0 0 3 -2             (34) and rank                         0 2 6π 2 0 1 0 1 2 0 2 0 -1 0 -1 0                         = 3 < rank                                     0 2 6π 2 0 1 0 1 2 0 2 0 -1 0 -1 0                         0 0 3 -2                                     = 4 (35)
Carter [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] Mixed iterative Algorithm PRDV algorithm [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] ν int (rad) Figure 11 shows the in-plane trajectory of the chaser for Carter's solution (green) and mixed iterative algorithm (blue).

Interestingly, the simulation of Carter's maneuver planning as proposed in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] leads to an error at the final point of the rendezvous even with the Hill-Clohessy-Wiltshire state-transition matrix. The final academic case of Carter is recalled in the Table 5. As will be seen in the sequel, this simple example is particularly interesting since it exhibits numerical hurdles to find the optimal solution that is quite surprising for so simple an example.

Eccentricity e = 0 ν 0 0 rad X T 0 [ 1 0 0 0 ] ν f 2π rad X T f [ 0 0 0 0 ]

Analysis of Carter's solution

The solution proposed in [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] is a two-impulse solution with an initial and final thrust in opposite directions along the x axis of the LVLH frame. This solution is analytical in the sense that it is the minimum-fuel solution of the linear system:

R hcw (0)∆V (0) + R hcw (2π)∆V (2π) = φ -1 hcw (2π)X f -φ -1 hcw (0)X 0 = u f ( 36
)
and may readily be computed as:

∆V (0) = 1 6π 0 T ∆V (2π) = - 1 6π 0 T ( 37
)
This conjecture appears to be reasonable if the following result from [START_REF] Carter | Quadratic-based computation of four-impulse optimal rendezvous near circular orbit[END_REF] is recalled.

Theorem 4.1 ([18])

For HCW rendezvous problems, there is no optimal four-impulse conditions for boundary conditions defined as X T 0 =

x 0 z 0 0 0 , X T f = 0 0 0 0 .

Unfortunately, the solution (37) is suboptimal as may be demonstrated by the associated primer vector that does not satisfy all the necessary conditions of optimality. λ v (ν) is defined by:

λ v (ν) = R T (ν)λ ( 38
)
where λ is a real solution of the polynomial system (39):

R(0)R(0) T ∆ v 0 λ + R(2π)R(2π) T ∆ v f λ = u f λ T R(0)R(0) T λ = 1 λ T R(2π)R(2π) T λ = 1 (39)
minimizing λ T u f . This minimum-fuel solution λ * is given by λ * = -0.106103295393226 0.212206590785999 0.000443462832762 -0.99912

T ( 40
)
The norm and the locus of the primer vector associated to Carter's solution λ v (ν) c = R T (ν)λ * are depicted on Figure 12. We get that it violates the necessary condition of optimality (9) since its maximum is numerically given by:

λ c vm = max ν∈[ν 1 ,ν f ] λ v (ν) c = 1.052327525631438 (41)
The sub-optimality of Carter's solution may be formally proved by considering that λ * must be the solution of the following linear system:

            36π 2 + 8 4 12π -6π 4 2 0 0 12π 0 8 -4 -6π 0 -4 2                         λ 1 λ 2 λ 3 λ 4             =             -6π 0 0 0             (42)
Hence, λ * may be parametrized as:

λ * = - 1 3π 2 3π λ 3 -1 + 2λ 3 T (43)
If λ * is a solution verifying the optimality condition (9) then

λ v (ν) c = R T hcw (ν)λ * < 1, ∀ ν ∈ (0, 2π) (44)
The condition (44) is equivalent to the existence of λ 3 ∈ R such that ∀ ν ∈ (0, 2π):

λ v (ν) 2 = R T hcw (ν)λ 2 < 1 = 4(cos ν -1) 2 + sin 2 ν λ 2 3 - 4 π (cos ν -1) [ν -π -sin ν] λ 3 + (3(ν -π) -4 sin ν) 2 + 4(cos ν -1) 2 9π 2 < 1 (45)
For ν = π 20 , it is easy to show that the second order polynomial of (45) is always strictly greater than 1, ∀ λ 3 ∈ R. Carter's solution is therefore not optimal.

Analytical two-impulse solution

Let us now consider this problem but without setting the two times of thrusting a priori but by defining (ν 1 , ν 2 ) as decision variables. The two-point boundary value problem reads as:

X f -Φ(2π, 0)X 0 = Φ(2π, ν 1 )B∆V (ν 1 ) + Φ(2π, ν 2 )B∆V (ν 2 ) (46)
and is equivalent to the system (47) to be solved:

            -6π + 3ν 1 -4 sin ν 1 2 -2 cos ν 1 -6π + 3ν 2 -4 sin ν 2 2 -2 cos ν 2 -2 + 2 cos ν 1 -sin ν 1 -2 + 2 cos ν 2 -sin ν 2 -3 + 4 cos ν 1 -2 sin ν 1 -3 + 4 cos ν 2 -2 sin ν 2 2 sin ν 1 cos ν 1 2 sin ν 2 cos ν 2                         ∆V x (ν 1 )
∆V z (ν 1 )

∆V x (ν 2 ) ∆V z (ν 2 )             =             -1 0 0 0             M(ν 1 , ν 2 )∆V = ϖ f (47) The determinant of M(ν 1 , ν 2 ) is: det M(ν 1 , ν 2 ) = 16 sin 2 ((ν 2 -ν 1 )/2) -3(ν 2 -ν 1 ) sin(ν 2 -ν 1 ) = g(θ ) (48)
where 0 ≤ θ = ν 2 -ν 1 ≤ 2π. The function g(θ ) is a strictly positive on (0, 2π) and vanishes for 0 and 2π as may be easily verified [START_REF] Prussing | Optimal multiple-impulse orbital rendezvous[END_REF]. Cases θ = 0 and θ = 2π may be excluded for obvious reasons. We get the vector of thrusts as functions of (ν 1 , ν 2 ).

            ∆V x (ν 1 ) ∆V z (ν 1 ) ∆V x (ν 2 ) ∆V z (ν 2 )             = M -1 (θ )ϖ f = 1 g(θ )             -sin θ -4 sin 2 (θ /2) sin θ -4 sin 2 (θ /2)             (49)
The L 1 /l 2 minimum-fuel rendezvous problem is nothing but the parametric optimization problem (50) where θ is the only decision variable. min θ 2 sin 2 θ + 16 sin 4 (θ /2)

16 sin 2 (θ /2) -3θ sin θ = J l 2 (θ ) s.t. 0 < θ < 2π (50)
As it is, this problem appears difficult to solve analytically. In fact, (50) is a convex optimization problem for which a minimum may be computed via the computation of the only zero of its derivative (51) on the interval (0, 2π) by a Newton method.

dJ l 2 (θ ) dθ = -2 18θ + 8 sin(2θ ) -9 sin(3θ )/4 -37 sin θ /4 + 24θ (2 sin 2 (θ /2) -1) -6θ (2 sin 2 (θ ) -1) (16 sin 2 (θ /2) -3θ sin θ ) 2 16 sin 4 (θ /2) + sin 2 θ (51)
We get the minimizer θ * of (50) as

θ * = 6.230033575529312 (52)
The first and second derivatives of g(θ ) are shown on Figures 13a and13b confirming the convexity of the problem on the interval (0, 2π).

d 2 J l 2 (θ ) dθ 2
Fig. 13:

dJ l 2 (θ ) dθ and d 2 J l 2 (θ ) dθ 2 on (0, 2π)
Finally, the optimal cost and the optimal maneuvers for Carter's example are given by:

J l 2 (θ * ) = 0.105954087364712             ∆V x (ν 1 ) ∆V z (ν 1 ) ∆V x (ν 2 ) ∆V z (ν 2 )             =             -0.052902333870518 -0.002812512822111 0.052902333870518 -0.002812512822111             (53)
These results show that the optimal two-impulse solution for this example is non unique. There is an infinity of optimal solutions with the same consumption and a couple of optimal impulse times verifying

(ν * 1 , ν * 2 ) ∈ (0, 2π -θ * )×(θ * , 2π) with ν * 2 -ν * 1 = θ * .

The mixed iterative algorithm solution

The results of the mixed iterative algorithm are presented in the Table 6 and clearly show that Carter's solution is indeed not an optimal solution. The mixed algorithm converges after 13 iterations to a four-impulse solution illustrated by the Figure 15a.

Mixed iterative Algorithm ν i (rad) 0 0.0681 6.1837 2π

∆V (ν 0 ) T [ 0.021453 0.001138 ] ∆V (ν 1 ) T [ 0.031449 0.001685 ] ∆V (ν 2 ) T [ 0.006797 0.000364 ] ∆V (ν f ) T [ 0.046105 0.002446 ] L 1 /l 2 Cost 0.105954
Table 6: Results of the mixed iterative algorithm for Carter's first example [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] The norm of the primer vector is given on Figure 14 where λ vm (ν) = 1.0000048. The distribution of this solution may be defined as two pairs of almost simultaneous impulses at the beginning and at the end of the rendezvous (as is illustrated on the plot of the trajectory in the orbital plane at Figure 15b).

To save fuel and allow for in-flight testing throughout the FFIORD experiment, the rendezvous maneuver must last several orbits. Duration of the rendezvous is approximately 12 hours for an expected average cost of 20 cm/s [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF].

The mixed iterative algorithm achieves optimization within 13 seconds and within 3 iterations. Global optimality of this three-impulse solution has been established in [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] [START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF] for the PRISMA case study Figure 16b shows primer vector magnitude during transfer. Note the low magnitude of the second impulse (0.002 m/s) with respect to the initial and final velocity increments (0.0492 m/s and 0.0513 m/s) but these velocity increments play a significant role in the optimality of the result. In particular, they provide the right chaser orientation for the long drift (61400 s) between the second impulse and the final one. Indeed, the designer could be tempted to remove this interior impulse and resort to the suboptimal two-impulse strategy. The latter solution proves to be strongly suboptimal since its L 1 /l 2 cost is 27 % greater than the optimal solution (0.14506 m/s). The long drifting period of 61400 s of the optimal solution is clearly illustrated in Figure 16a where the in-plane trajectory and impulse positions are represented. Finally, it is worth noticing that the optimal cost is half the expected average cost of 20 cm/s [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF]. A new numerical algorithm based on heuristic rules deduced from the work of [START_REF] Lion | Primer vector on fixed-time impulsive trajectories[END_REF] and tools from algebraic geometry has been proposed to address the issue of time-fixed optimal rendezvous in a linear setting. This algorithm is a mixed iterative algorithm optimizing over the number of impulses with a low numerical complexity mainly consisting in the solution of a small size polynomial system of equations. Even without a formal proof of convergence, this heuristic algorithm appears to be very efficient in practice on very different missions ranging from circular to elliptic rendezvous missions. In particular, it is not necessary to resort to some local optimization scheme as in the references [START_REF] Jezewski | An efficient method for calculationg optimal free-space n-impulse trajectories[END_REF] and [START_REF] Jezewski | Primer vector theory applied to the linear relative-motion equations[END_REF] to eliminate cusp occurrences in the graph of the norm of the primer vector.

Despite the good numerical results presented, some improvement can still be expected if more sophisticated transition matrices including orbital perturbation effects are used such as atmospheric drag or gravity harmonics term J 2 .Another avenue of research deals with the extension of this algorithm for optimal trajectory planning with collision avoidance constraints.

| ν ν 1 ν f ν m 1 (Fig. 1 :

 11 Fig. 1: Nonoptimal primer vector norms

  Iteration 3: ν b = ν m = 1.7331

Fig. 5 :

 5 Fig.5: Details of iterations 4-6 for the third case of Carter[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]

  Figure 6b depicts the trajectory of the chaser in the orbital plane. Only linear simulations based on Hill-Clohessy-Wiltshire transition matrix have been performed.

Fig. 6 :

 6 Fig. 6: State and control trajectories for the Carter's third example

1 : 3 : 8 Fig. 8 : 5 :Fig. 9 : 7 :Fig. 10 :

 138859710 Fig.8: Iterations 1-3 details for Carter's case study 2[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]

Fig. 11 :

 11 Fig. 11: In-plane trajectory and impulse dates for Carter's case study 2: Mixed iterative algorithm (blue) and Carter's solution (green)

Fig. 12 :

 12 Fig. 12: Norm and primer vector locus for Carter's two-impulse solution

Fig. 14 :

 14 Fig. 14: Norm of the primer vector and cost evolution

Fig. 16 :

 16 Fig. 16: Primer vector norm and trajectory for the PRISMA mission

Table 2 :

 2 Results comparison for the third case of[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]

	]

Table 4 :

 4 Results from[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF], mixed iterative and PRDV algorithms for Carter's case study 2[START_REF] Arzelier | Using polynomial optimization to solve the fuel-optimal impulsive rendezvous problem[END_REF]

	∆V (ν 0 ) T ∆V (ν 1 ) T ∆V (ν f ) T	π 2 ≃ 1.5708 [ 1.6294 -0.6667 ] [ 0.3901 0.0964 ] [ -0.0633 -0.0259 ]	2.4119 [ 1.7775 -0.3828 ] [ 0.2896 -0.0165 ] [ -0.0672 -0.0143 ]	2.4085 [ 1.7771 -0.38449 ] [ 0.28995 -0.015971 ] [ -0.06706 -0.014384 ]
	L 1 /l 2 cost	2.2307	2.1770	2.1772

Table 5 :

 5 Data for Carter's first example[START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF]

Table 8 :

 8 by running the PDRV algorithm. Results of the mixed iterative algorithm and PRDV algorithm

		PRDV Algorithm	Mixed iterative algorithm
	t int (s) ν int (rad) ∆V (ν 0 ) T ∆V (ν 1 ) T ∆V (ν f ) T	3189.3 3.4285 [ -0.04911 0.002152 ] [ -0.002038 0.0000099 ] [ -0.002039 0.0000112 ] 3198.6 3.4377 [ -0.04911 0.001933 ] [ 0.051315 0.001423 ] [ 0.051316 0.001404 ]
	L 1 /l 2 Cost m/s	0.102525	0.10252

This work was supported by CNES Grant R-S07/VF-0001-065.

Fig. 15: State and control trajectories for the Carter's first case [START_REF] Carter | Optimal impulsive space trajectories based on linear equations[END_REF] Keeping in mind that the optimal solution has been computed analytically as a two-impulse solution with one initial and one final coasting period and recalling the result presented by theorem 4.1, the result of the mixed iterative algorithm may be analyzed as a tight approximation of the genuine optimal two-impulse solution.

Case study 4

Following the first three academic numerical examples, a more realistic illustration based on PRISMA [START_REF] Berges | CNES approaching guidance experiment within FFIORD[END_REF] is now presented.

PRISMA programme is a cooperative effort between the Swedish National Space Board (SNSB), the French Centre National d'Etudes Spatiales (CNES), the German Deutsche Zentrum für Luft-und Raumfahrt (DLR) and the Danish Danmarks Tekniske Universitet (DTU) [START_REF] Larsson | Fuel Efficient Relative Orbit Control Strategies for Formation Flying and Rendezvous within PRISMA[END_REF]. Launched on June 15, 2010 Yasny (Russia), it was intended to test in-orbit new guidance schemes (particularly autonomous orbit control) for formation flying and rendezvous technologies. This mission includes the FFIORD experiment led by CNES which features a rendezvous maneuver (formation acquisition). The orbital elements of the target orbit as well as initial and final rendezvous conditions are listed in Table 7. Table 7: PRISMA rendezvous characteristics


J J Stirnemann

F Comte

DENSITY ESTIMATION OF A BIOMEDICAL VARIABLE SUBJECT TO MEASUREMENT ERROR USING AN AUXILIARY SET OF REPLICATE OBSERVATIONS

Keywords: Density estimation, non-parametric methods, measurement error, gestational age, replicate measurements

Correcting for measurement error the density of a routinely collected biomedical variable is an important issue when describing reference values for both healthy and pathological states. The present work addresses the problem of estimating the density of a biomedical variable observed with measurement error without any a priori knowledge on the error density. Assuming the availability of a sample of replicate observations, either internal or external, which is generally easily obtained in clinical settings, an estimator is proposed based on non-parametric deconvolution theory with an adaptive procedure for cut-off selection, the replicates being used for an estimation of the error density. This approach is illustrated in two applicative examples: i) the systolic blood pressure distribution density using the Framingham Study dataset and ii) the distribution of the timing of onset of pregnancy within the female cycle, using ultrasound measurements in the first trimester of pregnancy.

Introduction

The variability of a biomedical variable in a target population depends on several factors. However, its measurement inevitably exposes to measurement error, which adds to the observed variability. This is particularly important in clinical settings or large scale epidemiological studies, when dealing with routinely collected variables, like blood pressure or kaliema, the measurement of which relies on simple devices and/or depends on the circumstances and the operator. Depending on the magnitude of the measurement error, a direct consequence is that the probability distribution function of the observable variable deviates from the probability distribution function of the unobservable true underlying variable. Therefore, measurement error will alter any statistical inference depending on the amount of error the true measurement is contaminated with. In epidemiology and regression for example, measurement error in the exposure variable will bias the result towards the null hypothesis. Numerous examples of methods for dealing with such an error have been suggested in epidemiology (see [START_REF] Budtz-Jørgensen | Consequences of exposure measurement error for confounder identification in environmental epidemiology[END_REF] for an example in environmental epidemiology and [START_REF] L S Freedman | A comparison of regression calibration, moment reconstruction and imputation for adjusting for covariate measurement error in regression[END_REF] for an example in a dietary intake cohort study). Similar considerations have been discussed in clinical trials, mostly when the design involves surrogate end-points and biomarkers [START_REF] Sarkar | Quantifying the treatment effect explained by markers in the presence of measurement error[END_REF], [START_REF] Li | Adjustment for the measurement error in evaluating biomarkers[END_REF]. Other common examples include genetics and population association studies, involving both clinical [START_REF] Lobach | Haplotype-based regression analysis and inference of case-control studies with unphased genotypes and measurement errors in environmental exposures[END_REF], [START_REF] Barendse | The effect of measurement error of phenotypes on genome wide association studies[END_REF] and technical [START_REF] Bergemann | Signal quality measurements for cDNA microarray data[END_REF] aspects of variables with measurement error.

In the present work, we are interested in the basic problem of estimating the density of a variable when it is measured with an unknown error. There are numerous examples that justify the need for an accurate estimation of the density of a biomedical variable such as the establishment of reference ranges. More particularly, taking potential measurement noise into account should improve the threshold values when planning a screening program based on blood biomarkers. In descriptive epidemiology, removing the measurement error might help in uncovering a bimodal distribution of the variable of interest, thus suggesting a mixture of two different subpopulations. However, since the error is not directly observed, the issue of estimating the density of the error-free unobserved variable is not straightforward. Depending on the setting, several solutions exist for estimating density of the true variable: i) when the density of the error is assumed to be known, classical deconvolution algorithm may be used to achieve a non-parametric estimation using kernel [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], [START_REF] Liu | A consistent nonparametric density estimator for the deconvolution problem[END_REF], [START_REF] Stefanski | Deconvoluting kernel density estimators[END_REF], [START_REF] Hesse | Data-driven deconvolution[END_REF], [START_REF] Delaigle | Bootstrap bandwidth selection in kernel density estimation from a contaminated sample[END_REF] or wavelet methods [START_REF] Fan | Wavelet deconvolution. Information Theory[END_REF], [START_REF] Pensky | Adaptive wavelet estimator for nonparametric density deconvolution[END_REF]. With biomedical variables, this may occur when gross inspection of the data strongly suggests a Gaussian error. However, in many circumstances Gaussian approximation of the error may not hold, as exemplified later. Furthermore, when the parameters of the error distribution are unknown, these methods require a sample estimation of both the mean and variance of the error as an additional estimation step for plugging into a density estimation algorithm. Many of these efficient algorithms are presently available in the R statistical decon package [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF], [START_REF] Wang | Deconvolution estimation in measurement error models: The R package decon[END_REF]; ii) the density of the error may be unknown but a sample of pure errors without signal is available. In this circumstance, deconvolution algorithms have been implemented and have shown good properties regarding convergence Diggle and [START_REF] Hall | A fourier approach to nonparametric deconvolution of a density estimate[END_REF], [START_REF] Neumann | On the effect of estimating the error density in nonparametric deconvolution[END_REF], [START_REF] Comte | Data-âĂŘdriven density estimation in the presence of additive noise with unknown distribution[END_REF], [START_REF] Wang | The effects of error magnitude and bandwidth selection for deconvolution with unknown error distribution[END_REF]. However, this situation is less likely to occur in biomedical research especially epidemiological and clinical data, since a sample of signal-free errors requires a specific experimental design, in which signal-free samples have some meaning and are available. However, this situation occurs when dealing with measurement devices such as biological assays or imaging intensity measures. Indeed, assessing the intrinsic characteristics of any new measuring technique is a mandatory step although not always fully informative for the real-life setting.

The present work focuses on a situation where no hypothesis regarding the measurement error distribution can be made and where there is no access to a sample of signal-free measurement errors. However, as we will see in the following Section, errors are considered random, additive and homoscedastic. An estimator of the density of the measurement error-free variable is proposed based upon a second sample comprising replicate observations obtained in the same way as the primary sample of single observations. This estimator has shown good asymptotic properties both theoretically and in simulations Comte et al. [2011]. To illustrate its potential use in biomedical data, we apply this method in two real-data case studies: i. Systolic blood pressure in the Framingham Study on coronary heart disease. Replicate data is used to estimate the distribution of systolic blood pressure, measured twice on the same visit. The estimate is compared to both a naive estimator and to an estimation assuming a Gaussian approximation of the error density. ii. Timing of onset of pregnancy within the female cycle in a cohort of pregnant women.

We are interested here in estimating the density of the time interval between last menstrual period and onset of pregnancy. Since the true date of onset of pregnancy is unknown, we may only observe a noisy observation of this time interval. A first trimester ultrasound measurement of crown-rump length of the embryo is used as a noisy observation for dating pregnancy. Replicate observations are obtained from a sample of twin pregnancies using ultrasound measurements on each twin.

In Section 2 we introduce the general frame and notations of deconvolution models followed by our estimator with replicate measurements. Case-studies are presented in Sections 3 and 4. We finish with some discussion in Section 5.

Deconvolution with replicate measurements

Only the outline of the model and methods is presented here. Technical details regarding theoretical developments and convergence properties may be found in [START_REF] Comte | Deconvolution estimation of onset of pregnancy with replicate observations[END_REF] Comte et al. [2011].

2.1. Model and notations. We will denote Y j the noisy observation and X j the unobserved variable for subject j. The goal is to estimate f the density of X. The classical formulation of measurement-error models yields:

(1)

Y j = X j + ε j , j = 1, . . . , n
where ε j is an error term identically and independently distributed with an unknown density f ε . Furthermore, we consider the specific case of a homoscedastic error: the sequences (ε j ) 1≤j≤n and (X j ) 1≤j≤n are assumed independent. The density of Y j denoted f Y is the convolution of the densities of X j and ε. Denoting the convolution operator, we have:

(2)

f Y (x) = (f f ε )(x) = f (x -u)f ε (u)du
Equivalently, taking the characteristic function of each density denoted by an asterisk, we have:

(3)

f * Y (u) = f * (u) × f * ε (u) Since Y is observed, a natural estimator of f *
Y will be the empirical characteristic function. However, in most biomedical circumstances, the density of the error will be unknown prior to analysis and a sample of pure errors will not be available and sometimes altogether nonrealistic. On the contrary, it is often possible to obtain a sample of repeated measurements, which we call replicates throughout this article not to confuse the reader with longitudinal observations. We define replicates as two or more error-contaminated observations of same nature (with identically distributed errors) of a single unobserved true quantity. This is often the case in reproducibility studies for example, when two or more observers will measure the same biological parameter. We suggest that an estimator of the density of an error-contaminated variable can be obtained using replicate observations in this circumstance.

Estimators and deconvolution with replicate measurements.

Let us assume we have a sample of size M of 2L replicate noisy observations of ε k,2 -1 andε k,2 , for k = 1, . . . , M , = 1, . . . , L, independent and identically distributed. We assume the sequences ( X k ) 1≤k≤M , (ε k,2 -1 ) 1≤k≤M and(ε k,2 ) 1≤k≤M are independent. Therefore, we consider that two independent samples are available: the first, of size M , containing replicate observations and the second, of size n containing nonreplicate observations. Density estimation by deconvolution with replicate observations has been previously studied [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF], [START_REF] Li | Nonparametric estimation of the measurement error model using multiple indicators[END_REF], [START_REF] Meister | Deconvolution from non-standard error densities under replicated measurements[END_REF], showing that an estimator may be achieved using the sample of replicates to estimate the noise density and to improve the deconvolution step. Although related to the estimator given by [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF], we suggest a new estimator that has shown good finite sample and asymptotic properties Comte et al. [2011]. Aside from the assumption of independence, a reasonable assumption is that ε is symmetric and that its characteristic function never vanishes. This assumption implies that the characteristic function is real-valued and strictly positive. Under this hypothesis, observing that

X k : (4) Y k,2 -1 = X k + ε k,2 -1 , Y k,2 = X k + ε k,2 , k = 1, . . . , M, = 1, . . . , L with X k ,
Y k,2 -1 - Y k,2 reduces to the difference of two independent, identically distributed random variables ε k,2 -1 -ε k,2 , we have the following relationship between the corresponding characteristic functions f * Y k,2 -1 -Y k,2 (u) = (f * ε (u)) 2 Under the assumption of a symmetric error, f * ε is real-valued.
Hence an estimator of the square of f * ε is obtained by taking the real part of the estimator of

f * Y k,2 -1 -Y k,2 : (5) (f * ε ) 2 (u) = 1 M L M k=1 L =1 cos u(Y k,2 -1 -Y k,2 )
We also define an estimator for f * Y which is the empirical characteristic function of the independent noisy observations (Y j ) 1≤j≤n and (Y k,1 ) 1≤k≤M :

(6) f * Y (u) = 1 n + M   n j=1 e iuY j + M k=1 e iuY k,1  
Because of numerical tractability, (5) cannot be plugged as such in (3) together with (6). In order to prevent the effect of dividing by small numbers in (3), we define a truncated estimator of f * ε , as suggested by [START_REF] Neumann | On the effect of estimating the error density in nonparametric deconvolution[END_REF][START_REF] Neumann | On the effect of estimating the error density in nonparametric deconvolution[END_REF]:

(7) 1 f * ε (u) = 1{ (f * ε ) 2 (u) ≥ (LM ) -1/2 } (f * ε ) 2 (u)
Finally, plugging ( 7) and ( 6) in (3) and using Fourrier inversion with a πm cut-off, we have the final estimator of f using a πm cut-off for integrability purpose :

(8) fm (x) = 1 2π πm -πm e -ixu f * Y (u) f * ε (u) du
This estimator can also be viewed as a deconvolution kernel estimator with bandwidth 1/(πm). Comte et al. [2011]. In brief, the estimator is projected on an orthonormal basis built from sinus cardinal functions defined by ϕ(x) = sin(πx)/πx. Such an orthonormal basis allows a smooth estimate while offering desirable computational properties. Concurrently, the choice of the cut-off πm is based upon an estimation of an optimal value for m defined by a bias-variance compromise:

m opt = argmin m f -f m 2 + var(m)
where f -f m 2 is the bias between the true unknown density function f and the function that is estimated by (8) while var(m) is the unknown variance of the estimator, function of m. Each of these two quantities may be estimated and plugged in to estimate an appropriate cut-off for the estimator.

The final algorithm allows a fully automatized estimation of f with cut-off selection and has been implemented in R. Compared to the estimator presented in [START_REF] Delaigle | On deconvolution with repeated measurements[END_REF], the presented method provides an adaptative selection procedure for the optimal bandwidth πm and uses a different orthonormal basis for projection. As presented in Comte et al. [2011], comparative simulations show that the presented estimator outperforms the earlier in terms of pointwise risk.

Case study: distribution of systolic blood pressure in the Framingham data

To exemplify the method in a simple case study, we consider data from the Framingham Study on coronary heart disease discussed in [START_REF] Carroll | Measurement Error in Nonlinear Models: A Modern Perspective, Second Edition[END_REF] [START_REF] Carroll | Measurement Error in Nonlinear Models: A Modern Perspective, Second Edition[END_REF] and in [START_REF] Wang | Deconvolution estimation in measurement error models: The R package decon[END_REF] [START_REF] Wang | Deconvolution estimation in measurement error models: The R package decon[END_REF] where systolic blood pressure (SBP) were measured. An accurate denoised estimation of SBP may be of interest in several ways: population reference ranges could be substantially refined; the results of epidemiologic studies or trials using systolic blood pressure as a covariate or as an outcome could be improved; considering the measuring devices are unbiased, the denoised estimate is independant of the measuring technique and may be regarded as a gold-standard since any device-specific noise has been removed.

The data consist of measurements of SBP in 1,615 males on several visits on an 8-year follow-up. During each visit, each patient had its blood pressure measured twice. Since each measurement is contaminated with an unknown error due to either human, biological or technical variations, the density of blood pressure is unknown. Considering only the first visit, our goal is to estimate the density of SBP , using the replicate measurements denoted SBP 1 and SBP 2 on each individual. In these data, a plot of SBP 1 -SBP 2 shows that the hypothesis of a Gaussian error would roughly hold in this case as shown in Figure 1 and standard packaged deconvolution algorithms could be used. However, for the sake of the experience, let's not assume a parametric known distribution for the error. To fit our setting, we split the sample into two separate datasets: one containing the first 500 replicate observations of SBP 1 and SBP 2 and the second containing only the 1115 last observations of SBP 2. The hypothesis of homoscedastic error is investigated graphically by plotting SBP 1 -SBP 2 versus (SBP 1 + SBP 2)/2 in Figure 1 showing no significant relationship. Deconvolving SBP 2 using the replicate observations yields the estimate presented in Figure 2 (solid line). This estimate is compared to the deconvolution estimator obtained when considering the noise as Gaussian (dotted line) using the decon package in R Wang and [START_REF] Wang | Deconvolution estimation in measurement error models: The R package decon[END_REF]. Despite ungraceful ripples in the tails of our estimate, our estimate (solid line) shows a more peaked mode compared to the estimate with Gaussian error. This finding demonstrates that the hypothesis of a Gaussian error may in fact hamper the estimation since it seems to underestimate the amount of noise present in the data. Furthermore, the distributional mode is slightly right-shifted using our estimator. Compared to both deconvolution estimators, the naive kernel estimator (dashed-line) of the raw SBP 2 data underestimates the peak since the unknown error adds variance to the true distribution.

Case study: estimation of density of onset of pregnancy

Except in the specific case of in vitro fertilization, the precise date of pregnancy is unknown in women conceiving naturally. Therefore, it must be estimated using proxy measures that are subject to measurement error (see [START_REF] D B Dunson | Day-specific probabilities of clinical pregnancy based on two studies with imperfect measures of ovulation[END_REF][START_REF] D B Dunson | Day-specific probabilities of clinical pregnancy based on two studies with imperfect measures of ovulation[END_REF] and [START_REF] D B Dunson | Assessing human fertility using several markers of ovulation[END_REF][START_REF] D B Dunson | Assessing human fertility using several markers of ovulation[END_REF] for a comprehensive study regarding measurement error in markers of ovulation). Ultrasound biometry of the embryo in early pregnancy has been proven the most accurate method for dating pregnancy in clinical practice Gardosi et al. [1997b,a], [START_REF] Mongelli | Birth weight, prematurity and accuracy of gestational age[END_REF], [START_REF] Tunón | A comparison between ultrasound and a reliable last menstrual period as predictors of the day of delivery in 15,000 examinations[END_REF]. Numerous regression formulas have related early ultrasound biometry and gestational age both in spontaneously conceived pregnancies and pregnancies conceived with assisted reproduction techniques [START_REF] Sladkevicius | Ultrasound dating at 12-14 weeks of gestation. a prospective cross-validation of established dating formulae in in-vitro fertilized pregnancies[END_REF]. The most widely used formula is that of [START_REF] Robinson | Sonar measurement of fetal crown-rump length as means of assessing maturity in first trimester of pregnancy[END_REF] using crown-rump length in the first trimester of pregnancy. Therefore, for women conceiving spontaneously, early pregnancy ultrasound biometry provides a relevant noisy observation Y j of X j .

4.1. The data. The first dataset is a sample of singleton pregnancies with one noisy observation Y j of X j ; the second is a sample of replicate observations obtained from twin pregnancies yielding noisy observations Y k1 and Y k2 of X k .

4.1.1. Single observation sample. The data comprises basic information recorded at time of first trimester ultrasound and is provided by a single early pregnancy screening center in Paris, France. Since screening is part of a national policy, women are unselected and are likely to be a representative sample of the general population of pregnant women. First trimester screening visits are planned at around 12 weeks following last menstrual period (about 10 weeks of gestation) and between 11 and 14 weeks as determined by ultrasound.

Recorded data includes maternal age, characteristics of the cycles, last menstrual period and ultrasound biometric measurements in the fetus.

Cases were selected if conception was natural (i.e. without assisted reproductive techniques) and if the woman could recall last menstrual period. Because it could influence the recollection of last menstrual period, women were also excluded if ultrasound dating of the pregnancy had been performed prior to the visit.

Over a one-year period, 1,706 cases met these criteria. Menstrual cycles were considered regular in 1386 women whereas 320 women had irregular cycles. For each woman, the observed time interval between the last menstrual period and the date given by the ultrasound measurement of crown-rump length (CRL) was computed, representing the noisy observation Y j discussed in Section 2.1. 4.1.2. Replicate observation sample. Replicate independent observations are obtained from 86 spontaneous twin pregnancies in the same clinical setting as singleton pregnancies. Since onset of pregnancy is the same for both twins, the difference in CRL between the twins may be considered a measurement error. In more technical terms, we have shown in ( 5), that the empirical characteristic function of the difference of replicate noisy observations is an estimator of the square of the characteristic function of the noise itself. Therefore, each twin pregnancy has two observed time-intervals between LMP and DP based upon ultrasound, thus defining replicates of noisy observations Y k1 and Y k2 defined in Section 2.2.

It would be arguable here to check if the assumption of normality regarding ε j would hold, therefore allowing standard deconvolution algorithms with f ε considered as known. Denote σ2 the empirical estimator of the variance of Y k1 -Y k2 . Since both twins are exchangeable, we chose to compare the square of the standardized difference ((Y k1 -Y k2 )/σ) 2 to a χ 2 (df = 1) distribution. Graphically, Figure 3 shows that the hypothesis of a normal error is difficult to sustain, mostly because of departure in the quantile-quantile plot. Therefore, standard deconvolution methods will not apply since we cannot guess any further which potential density the error ε j originates from. Thus, the only way around is to consider the density of the error unknown as discussed in Section 2. 4.2. Estimation. We first consider only women with regular cycles. Deconvolution estimation of f in this population is presented in Figure 4. On the X-axis, t = 0 corresponds to the last menstrual period. The mode of the distribution is reached at 13 days meaning that in pregnant women with regular cycles, pregnancy is most likely to occur at that date within the cycle. Noteworthy, this distribution is skewed with a very low probability of pregnancy before 7 days and after 28 days. Therefore, even in women with regular cycles, there is a wide variation in the onset of pregnancy within a female cycle. Compared to the naive estimator not taking into account the error using standard kernel estimation routines (dashed line), our estimate shows a narrower and higher peak because of noise removal. 4.3. Influence of covariates. Several covariates have been shown to affect menstrual cycle characteristics including age, ethnicity, smoking, alcohol and other toxic substances

(Y k1 -Y k2 ) 2 /σ 2 versus the average (Y k1 -Y k2 )/2.
No significant trend was found using a nonparametric LOESS scatterplot smoother (solid line). [START_REF] Liu | Factors affecting menstrual cycle characteristics[END_REF]. The influences of age and cycle regularity are studied here as potential covariates for the timing of pregnancy. regular and irregular cycles is presented in the bottom panel of Figure 5. Consistent with Wilcox, Dunson and[START_REF] Wilcox | The timing of the "fertile window" in the menstrual cycle: day specific estimates from a prospective study[END_REF], the density of onset of pregnancy is shifted towards longer time intervals with a wider variation than women with regular cycles.

Concluding remarks

We have presented a general method for estimating the density of a variable with error in measurement. This method is conceived for frequent situations occurring in biological and clinical research settings and should allow researchers to improve analyses and provide more meaningful results. We emphasize that only reasonable hypotheses are required. In particular, estimation does not require a prior hypothesis regarding a parametric density of the error. Such parametric distributions for errors usually include Laplace and Gaussian densities which may often be an undesirable approximation. Owing to replicate observations, only the hypothesis of Gaussian errors may be easily investigated since the convolution of other distributions is unlikely to have a known parametric solution, except in very specific cases. The only hypothesis regarding distribution of the errors is symmetry, which seems a reasonable assumption in many situations involving biomedical data. We wish to emphasize that replicate observations are not independent: only the errors ε k1 and ε k2 must be independent and identically distributed as well as independent of the underlying true variable X k . Independence between errors in replicate observations is a classical and reasonable assumption. We acknowledge that independence between the error and the true variable, defining homoscedasticity, may be difficult to verify. However, following [START_REF] Bland | Measuring agreement in method comparison studies[END_REF] [START_REF] Bland | Measuring agreement in method comparison studies[END_REF], in most cases graphical inspection of a scatterplot of (Y k1 -Y k2 ) versus (Y k1 +Y k2 )/2 should be sufficiently conclusive. In case of evidence of heteroscedasticity, solutions have been recently suggested for replicate measurements under the assumption of a Gaussian error, considering the variances as random [START_REF] Mcintyre | Density estimation with replicate heteroscedastic measurements[END_REF].

Two case-studies were used to exemplify the use of deconvolution in analyses of biomedical data: in the Framingham data, we showed that using an unknown-error algorithm does not alter the estimation compared to classical algorithm assuming a fully known density for the error. In the pregnancy data, we conducted a more complete data analysis presenting novel results regarding timing of onset of pregnancy. The results of this second case study deserve further discussion regarding interpretation and potential applications. Using twin pregnancies may be arguable as a choice of replicate measurements since they may differ in some way from the population of singleton pregnancies. However, as far as we know, there does not seem to be a significant difference between singleton and twins for ultrasound dating of pregnancy [START_REF] Dias | Firsttrimester ultrasound dating of twin pregnancy: are singleton charts reliable?[END_REF]. Although it is related to previous results regarding the probability of ovulation during the menstrual cycle in healthy women often referred to as a "fertile window" [START_REF] D B Dunson | Day-specific probabilities of clinical pregnancy based on two studies with imperfect measures of ovulation[END_REF][START_REF] D B Dunson | Assessing human fertility using several markers of ovulation[END_REF], [START_REF] Wilcox | The timing of the "fertile window" in the menstrual cycle: day specific estimates from a prospective study[END_REF], the clinical question is different: we are interested in the timing of pregnancy rather than the timing of ovulation. These two events differ since pregnancy does not necessarily occur just following ovulation and because the potency of an ovulation to lead to a pregnancy may depend upon its timing itself. Therefore, the density of onset of pregnancy and the density of ovulation are likely to be different. Furthermore, dating ovulation can be performed accurately only in a specific prospective setting with intensive monitoring of biological samples throughout successive cycles, whereas we are concerned by the general population of pregnant women without enrollment in a specific experimental study. In a clinical perspective, the estimation of the density of timing of pregnancy is interesting in several ways: i) The variability of onset of pregnancy based upon last menstrual period displayed by the estimated density is much wider than the +/-5 days generally recognized as the interval of the error using ultrasound in the first trimester. Therefore, even in women with regular cycles and with a known date of last menstrual period, ultrasound will perform better than last menstrual period for dating pregnancy in clinical practice. ii) As discussed earlier, using any proxy for estimating onset of pregnancy in a single pregnant woman, such as crown-rump length ultrasound measurement, is subject to measurement error. However, in a Bayesian perspective and notwithstanding the larger variability in dating resulting from last menstrual period, the distribution of time to pregnancy may be used as a prior distribution for a likelihood function given by a proxy measurement. Thus, it may refine the estimate of a proxy measurement by incorporating last menstrual period as prior knowledge.

Figure 1 .

 1 Figure 1. Diagnostic plots for the Framingham data. Left-side panel: normal quantile-quantile plot of the standardized difference SBP 1 -SBP 2. Right-side panel: A scatterplot of the difference versus the average of the replicate observations for graphical inspection of potential heteroscedasticity. A nonparametric local scatterplot smoother (LOESS) is added (solid line).

Figure 2 .

 2 Figure 2. Estimation by deconvolution of the density of SBP 2 in the Framingham Study. The solid line represents the proposed estimator without assumption on measurement error, whereas the dashed line represents the deconvolution estimator with assumption of Gaussian measurement errors. The histogram of the raw observations of SBP 2 and the naive kernel estimator (dotted line) are given for comparison.

Figure 3 .

 3 Figure 3. Diagnostic plots in the pregnancy data. Left-side panel: quantile-quantile plot comparing (Y k1 -Y k2 ) 2 /σ 2 and a χ 2 (df = 1) distribution. Right-side panel: investigation of potential heteroscedasticity by plotting (Y k1 -Y k2 ) 2 /σ 2 versus the average (Y k1 -Y k2 )/2. No significant trend was found using a nonparametric LOESS scatterplot smoother (solid line).

Figure 4 .

 4 Figure4. Estimated density of onset of pregnancy for women with regular cycles (solid line). The histogram and naive estimator (dashed line) of the raw data Y j is given for comparison. The X-axis is the time interval between onset of pregnancy and last menstrual period (t=0).

  2.3. Outline of technical calculation of fm (x).Technical details regarding calculation of fm (x) can be found in[START_REF] Comte | Data-âĂŘdriven density estimation in the presence of additive noise with unknown distribution[END_REF], Comte, Rozenholc and Taupin[START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF] and[START_REF] Comte | Deconvolution estimation of onset of pregnancy with replicate observations[END_REF]. Details regarding the selection of an appropriate cut-off πm in (8) can be found inComte, Samson and Stirnemann

  4.3.1. Influence of maternal age. Maternal age was broken into four classes according to quartiles of increasing age in women with reported regular cycles. In each class of age, the density of onset of pregnancy was estimated independently. The results are presented in the top panel of Figure5. As expected, a noticeable shortening of the time interval between LMP and date of pregnancy is noticed with increasing age, consistent with the shortening of the first phase of menstrual cycle with age. In each class of age, the modes of the densities were 14.1 days, 13.1 days, 13 days and 13 days for women aged [18.1-29.3], ]29.3, 32.4], ]32.4,36] and ]36,48.1] respectively. The shape of the density also changes across age classes with positive skewness in younger women and negative skewness in older women consistent with the hypothesis of a shortening in the first phase of menstrual cycle with age. Noteworthy, the densities are more peaked in middle aged classes than in extreme classes showing that the probability of date of pregnancy is narrower around the mode in women between 29 and 36 compared to older and younger women. 4.3.2. Influence of regularity of cycles. Prior knowledge regarding individual menstrual characteristics may influence the dating of pregnancy since patients with irregular cycles are subject to more variation in timing of ovulation and therefore in timing of pregnancy as shown by[START_REF] Wilcox | The timing of the "fertile window" in the menstrual cycle: day specific estimates from a prospective study[END_REF] [START_REF] Wilcox | The timing of the "fertile window" in the menstrual cycle: day specific estimates from a prospective study[END_REF] in a study of fertility in non-pregnant fertile women. A comparison of timing of pregnancy in patients with reported

		0.12
		0.10
	Estimated density	0.06 0.08
		0.04
		0.02
		0.00
		0	5	10	15	20	25	30	35
					time (days)

Acknowledgements.

We thank Professor J.C. Thalabard for his motivating input as well as for his contribution in correcting the manuscript.