Roberto Armellin
email: armellin@aero.polimi.it

Pierluigi Di Lizia

Martin Berz

Kyoko Makino

Rigorous Global Optimization of Impulsive Planet-to-Planet Transfers in the Patched-Conics Approximation

Keywords: trajectory optimization, impulsive transfers, Taylor models, rigorous global optimization Trajectory optimization, impulsive transfers, Taylor models, rigorous global optimization

come   L'archive ouverte pluridisciplinaire

181x137mm (600 x 600 DPI)

Introduction

Interplanetary trajectory design problems are usually stated as optimization problems (Myatt et al. 2004, Di Lizia and[START_REF] Di Lizia | Advanced Global Optimisation for Mission Analysis and Design[END_REF] with the goal of maximizing the payload mass, or equivalently minimizing the fuel necessary on board. Interplanetary transfers are modeled either by continuously propelled arcs or by impulsive maneuvers. In the first case, e.g. when low-thrust propulsion systems are considered, the trajectory design is usually approached by applying optimal control theory, resulting in the solution of two point boundary value problems or large scale parametric optimizations [START_REF] Bryson | Applied Optimal Control[END_REF]

Taylor & Francis and I.T. Consultant [START_REF] Betts | Practical Methods For Optimal Control Using Nonlinear Programming[END_REF]. The second case is associated to chemical engine systems capable of high thrust magnitude but short firing times and low efficiency in terms of mass consumption. In this case, the thrust action is modeled as a sudden change in the spacecraft velocity and the optimization problem aims at finding the values of parameters describing the number and the magnitude of the velocity changes as well as the time instant in which they are applied. An impulsive transfer is generally simpler to be modeled and characterized by a reduced number of variables. Nevertheless, the solution of the optimization problem is challenging as the search space is usually huge and the objective function exhibits a large number of clustered minima. This feature even shows up in the simplest case of interplanetary transfer, namely a planet-to-planet transfer, and is even more emphasized when the gravitational pull of other celestial bodies -known as gravity assist maneuvers -and deep space maneuvers are considered. As local gradient-based methods converge to solutions corresponding to local minima, their applications to this category of problems is unsuitable. For this reason, a great effort has been recently spent to effectively prune unfavorable regions of the search space prior to the optimization run [START_REF] Izzo | Search space pruning and global optimisation of multiple gravity assist spacecraft trajectories[END_REF][START_REF] Bernelli-Zazzera | Global Trajectory Optimization: Can we Prune the Solution Space when Considering Deep Space Maneuvers[END_REF]) and to approach the problem from a global optimization point of view [START_REF] Yokoyama | Modified Genetic Algorithm for Constrained Trajectory Optimization[END_REF]Suzuki 2005, Vasile et al. 2005). The algorithms applied for this purpose are usually stochastic global optimizers of different kinds; differential evolution [START_REF] Storn | Differential Evolution-A Simple and Efficient Heuristic for global Optimization over Continuous Spaces[END_REF], genetic algorithms (Golberg 1989), and particle swarm optimizers [START_REF] Kennedy | Particle swarm optimisation[END_REF] being common examples (see [START_REF] Vinko | Benchmarking different global optimisation techniques for preliminary space trajectory design[END_REF] for a comparison between them). With these approaches, several runs of the same optimization code are necessary to assess the algorithm's efficacy, but there is no proof of gaining the global optimum of the problem at hand.

The application of a verified global optimizer to compute the verified global optimum of a simple planet-to-planet transfer is presented in this paper. By exploiting the capability of Taylor Model (TM) algebra [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF] to deliver the validated enclosure of functions, the optimizer returns the mathematically proven enclosure of the global optimum. As a consequence, the present work can be seen as a first step towards the verified global solution of complex impulsive interplanetary transfers.

The remainder of the paper is as follows. Section 2 is a brief presentation of Taylor models. Section 3 contains a description of the global optimizer COSY-GO. In Section 4 the validated solution of implicit equations is discussed in details. In Section 5 the mathematical formalization of the problem is described, focusing on the treatment of the ephemerides evaluation, Lambert's problem solution, and objective function evaluation within the TM framework. The discussion of the results achieved for Earth-Mars and Earth-Venus transfers follows, leaving conclusions and final comments to Section 6.

Notes on Taylor Models

Verified global optimization needs the determination of rigorous upper and lower bounds of the objective function in order to implement a branch and bound method [START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF]. The commonly used interval approach has excelled in solving this problem elegantly from both a formal and an implementational viewpoint. However, there are situations where the method has limitations for extended or complicated calculations because of the dependency problem, which is characterized by a cancellation of various sub-parts of the function that cannot be detected by direct use of interval methods. This effect often leads to pessimism and sometimes even drastic overestimation of range enclosure. Engineering Optimization 3

Furthermore, the sharpness of intervals resulting from calculations typically scales linearly with the sharpness of the initial discretization intervals. For complicated problems, and in particular higher dimensions, this sometimes significantly limits the sharpness of the resulting answer that can be obtained [START_REF] Makino | Efficient Control of the Dependency Problem Based on Taylor Model Methods[END_REF]).

The Taylor model approach enables the computation of fully mathematically rigorous range enclosures while largely avoiding many of the limitations of the conventional interval method [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF]. The method is based on the inductive local modeling of functional dependencies by a polynomial with a rigorous remainder bound, and as such represents a hybrid between formula manipulation, interval methods, and methods of computational differentiation (Berz et al. 1996, Griewank and[START_REF] Griewank | Automatic Differentiation of Algorithms: Theory, Implementation, and Application[END_REF].

An n-th order Taylor model of a multivariate function f that is (n + 1) times continuously partially differentiable on the domain D, consists of the n-th order multivariate Taylor polynomial P expanded around a point x 0 ∈ D representing a high order approximation of the function f , and a remainder error bound interval I for verification such that ∀x ∈ D, f (x) ∈ P (xx 0 ) + I.

(1)

Thus, the Taylor model (P, I) is a functional enclosure of f over D. From Taylor's theorem, it is clear that the width of the remainder interval I (i.e., the accuracy of the functional enclosure) can be chosen to scale with the domain size proportional to |xx 0 | n+1 . The practical computation of P and I is based on Taylor model arithmetic, which carries P and I through all the operations comprising I. (For a comprehensive textbook on interval analysis and rigorous computing the reader may refer to [START_REF] Moore | Introduction to Interval Analysis[END_REF]. By choosing the size |xx 0 | small and the order n sufficiently high, the size of the remainder interval I can be kept very small in practice. The bulk of the functional dependency is kept in the polynomial part P with point coefficients, and there is no interval arithmetic-associated inflation that happens in the polynomial part. Thus, the interval-related overestimation is rather optimally suppressed with the Taylor model method [START_REF] Makino | Efficient Control of the Dependency Problem Based on Taylor Model Methods[END_REF].

The implementation of the method in the code COSY Infinity (Makino 1998, Berz and[START_REF] Berz | COSY INFINITY Version 9 reference manual[END_REF]) supports binary operations and standard intrinsic functions, as well as the antiderivative operation which widens the applications of the method.

Altogether, the Taylor model approach has the following important properties:

(1) The ability to provide enclosures of any function given by a finite computer code list by a Taylor polynomial and a remainder bound with a sharpness that scales with order (n + 1) of the width of the domain.

(2) The computational expense increases only moderately with order, allowing the computation of sharp range enclosures even for complicated functional dependencies with significant dependency problem.

(3) The computational expense of higher dimensions increases only very moderately, significantly reducing the "curse of dimensionality".

The structure of Taylor models naturally represents a rich resource of information. In particular, the coefficients of the polynomial part P of a Taylor model are closely related to derivatives. That means when representing a function f by a Taylor model (P, I) on a computer, the local slope, Hessian and higher order derivatives are obtained almost free. When a task is focused on range bounding, those pieces of information become particularly useful.

While range bounding of Taylor Model with interval arithmetic in the naive sense Taylor & Francis and I.T. Consultant [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF] already exhibits superiority over mere interval arithmetic and the more advanced centered form [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF], the active utilization of those additional pieces of information in Taylor models has a lot of potential of developing efficient range bounders. Based on this observation, various kinds of Taylor-model-based range bounders have been developed [START_REF] Berz | Long-term stability of the Tevatron by verified global optimisation[END_REF], and among them the linear dominated bounder (LDB) and the quadratic fast bounder (QFB) are the backbones of Taylor model based verified global optimizer COSY-GO that will be discussed afterward.

The linear dominated bounder is based on the fact that, for Taylor models with sufficiently small remainder bound, the linear part of the Taylor model dominates the behavior; this is also the case for range bounding. The linear dominated bounder utilizes the linear part as a guideline for iterative domain reduction to bound Taylor models. Around an isolated interior minimizer, the Hessian of a function f is positive definite, so the purely quadratic part of a Taylor model (P, I) which locally represents f , has a positive definite Hessian matrix H. The quadratic fast bounder provides a lower bound of a Taylor model cheaply when the purely quadratic part is positive definite. More details on polynomial bounders are given by [START_REF] Makino | Verified Global Optimization with Taylor Model-based Range Bounders[END_REF].

COSY-GO

COSY-GO [START_REF] Berz | Long-term stability of the Tevatron by verified global optimisation[END_REF][START_REF] Makino | Performance of Taylor Model-based Algorithms for Unconstrained Global Optimzation[END_REF]) is a branch-and-bound optimization code employing local domain reduction techniques exploiting the bounding performances assured by TM methods. Should the global minimum of a sufficiently regular scalar function f on a given domain A ⊆ ℜ m wished to be evaluated, the algorithm starts with an initial value for the global optimum, the cutoff value, and then proceeds on analyzing at each step a subdomain for possible elimination or reduction. At each step the following tasks are performed (1) A rigorous lower bound l of the objective function is obtained on the subdomain of interest using various bounding schemes hierarchically with the hope of showing that l lies above the already established cutoff value, which will allow elimination of the subdomain. A first assessment is made whether the remainder bound of the Engineering Optimization ent is tested. For objective functions of nontrivial cost, as in the example at hand, also more sophisticated local searches within and near the current subdomain may be carried out.

The algorithm continues to reduce and examine the domain until the minimum dimension allowed is reached. The result of the optimization is the validated enclosure of the minimum of the problem.

Validated Solution of Implicit Equations

As will be shown in Section 5, the evaluation of the objective function for impulsive planet-to-planet transfers involves the solution of the implicit equations appearing in the ephemerides evaluation and Lambert's problem solution. The solution of these implicit equations must be carefully managed in the TM framework for a validated computation of the objective function. In this section the method for obtaining the validated solution of a scalar implicit equation is described first. Then, the extension to systems of equations is addressed, which is instrumental to the applications presented in Section 5.

Validated Solution of Scalar Implicit Equations

Let P (x) + I be an n-th order Taylor model of the (n + 1) times differentiable function

f over the domain D = [-1, 1] so that f (x) ∈ P (x) + I for all x ∈ D. (2)
Let R be an enclosure of the range of P (x) + I over D. Assume P has non-vanishing derivative everywhere in D. Without loss of generality, suppose P ′ (x) > d > 0 for all x ∈ D; i.e., according to the mean value theorem, P is strictly increasing in D.

Problem: Find a Taylor model c(y)

+ J on R such that any solution x ∈ D of the problem f (x) = y lies in c(y) + J. (3)
Without loss of generality, it is assumed that P (0) = 0. First, determine the polynomial c(y) = c 1 y + c 2 y 2 + . . . by inverting P . The inversion is performed by applying the algorithm illustrated by [START_REF] Berz | Modern Map Methods in Particle Beam Physics[END_REF], which reduces the inversion problem to the solution of an equivalent fixed point problem. Thus, the polynomial c(y) is the n-th order Taylor expansion of the inverse of f at the point 0 = P (0). For any given y, it returns an approximate value x = c(y) that satisfies P (x) ≈ y, and thus also f (x) ≈ y, depending on how well the inverse is representable by its Taylor expansion over the domain R. However, it is not the true inverse: evaluating P c(y) in the n-th order Taylor model arithmetic yields P c(y) ∈ y + J, where J is due to the terms of orders exceeding n in P c(y) ; and thus scales with at least order (n + 1).

Next a rigorous remainder J for c(y) is computed so that all solutions of f (x) = y lie in c(y) + J. To this aim the consequences of small corrections ∆x to c(y) are studied. According to the mean value theorem,

f c(y) + ∆x -y ∈ P c(y) + ∆x -y + I = P c(y) + ∆x • P ′ ξ(x) -y + I ⊂ y + J + ∆x • P ′ (ξ) -y + I = ∆x • P ′ (ξ) + I + J (4)
for some suitable ξ(x) that lies between c(y) and c(y) + ∆x. However, since P ′ is bounded below by d on [-1, 1], the set ∆x • P ′ (ξ) + I + J will never contain zero in [-1, 1], but outside of the interval

J = - I + J d , (5)
and thus no solution of f (c(y) + ∆x) = y exists there. Thus, the following theorem holds:

Theorem 4.1 Any solution x ∈ D = [-1, 1] of f (x) =
y lies in the Taylor model c(y) + J, with J given by (5).

Some remarks can be made.

1) No information on the derivative f ′ is needed, which is helpful in practice.

2) Indeed, if f is not monotonic, it is conceivable that c(y) + J contains multiple solutions of f (x) = y. 3) For any y ∈ [P (-1) + u(I), P (+1) + l(I)], at least one solution exists due to the intermediate value theorem. 4) For any y / ∈ [P (-1) + l(I), P (+1) + u(I)], no solution exists.

As an important consequence, the following corollary is obtained:

Corollary 4.2 If f is invertible over D and its Taylor polyomial has non-vanishing derivative over D, then the above c(y) + J is a Taylor model of order (n + 1) of its right inverse.

It is worth noting that, if f ′ (0) = 0, because of the inverse function theorem there is always a full neighborhood of 0 on which the requirements are satisfied.

The following example tests the above arguments on the sine function. The solutions of

sin(x) = y (6) for x ∈ D = [-1, 1] are of interest.
An exact representation of the solutions is available in analytical form, which is trivially given by

x = arcsin(y), (7)
where y lies in the range of sin(x) over D. Nevertheless, the above arguments are followed to find a Taylor model c(y)+J that encloses any solution of (6). As a consequence, c(y)+J must enclose the exact analytical solution (7). The Taylor model P (x) + I of the sine function on D is computed first. Then, P (x) is inverted to obtain c(y). Lastly, the interval J is computed to build the TM enclosure c(y) + J. The result of the procedure is shown in Figure 1 using 3-rd order expansions. The analytical solution is reported and compared with the TM enclosure. The range R of the Taylor model is assessed using the linear dominated bounder algorithm. As illustrated in the figure, for any y ∈ R, c(y) + J can be evaluated to identify a set that encloses the exact solution.

The accuracy of the enclosure depends on both the order of the Taylor expansions and the width of D. The effect of the order is investigated in Figure 2. Higher expansion orders yield sharper intervals J and, thus, higher accuracies. This is confirmed in figure, where the enclosures obtained using different orders are compared. It is also worth observing that the accuracy of the enclosure of the exact range of sin(x) over D, which is obtained using the LDB algorithm, improves with order.

The effect of the width of the domain D on the accuracy of the TM enclosure is analyzed in Figure 3. More specifically, the domain D is set as

D = 2 -k [-1, 1], (8)
and the width of J is computed and plotted for different k at fixed expansion orders. The logarithmic scale in the figure highlights a linear dependence with k, the slopes of the lines being directly related to n + 1. This confirms that the width of J scales with order n + 1 of the width of D.

Validated Solution of a System of Implicit Equations

The arguments introduced in Section 4.1 are extended in the following to the validated solution of a system of implicit equations. Consider the vector-valued function f (x) : (9)

D ⊂ R m → R m , with f i ∈ C n+1 (D),
Let R be an enclosure of the range of P (x)+I over D. Without loss of generality, assume

P (0) = 0 (10) J P (0) = I, (11)
where J P is the Jacobian of P and I is the identity matrix. Note that the assumption (10) can be met by adding a vector of constants to f , whereas (11) can be satisfied by an affine transformation in the codomain of f as long as the Jacobian of P is nonsingular at the origin. Let C(y) be the n-th order Taylor expansion of the inverse of P (x) about P (0). Similarly to the scalar case, C(y) is not the true inverse. Consequently, describing the range R with a TM number and evaluating P C(y) over R in n-th order Taylor model arithmetic yields

P C(y) ∈ y + J, ( 12
)
where J is due to the terms of orders exceeding n in P (C(y)), and thus scales with at least order (n + 1). Engineering Optimization 9 for all x ∈ D, with

Theorem 4.3 If ∂P i ∂x i (x) ≥ d > 0, ∀i ∂P i ∂x j (x) ≤ k, ∀j = i, (13)
d -(n -1) k > 0, (14)
then any solution x ∈ D of the problem f (x) = y lies in the Taylor model

C(y) + J, (15)
where J = [-a, a] n , with a defined as

a = max i=1,...,n |I i + Ji | d -k (n -1) . ( 16
)
The absolute value of an interval is defined as |I| = max |l(I)|, |u(I)| , where l(I) and u(I) are the interval lower and upper bounds.

Proof The theorem is proved by studying

f C(y) + ∆x -y, with ∆x such that C(y) + ∆x ∈ D. ( 17
)
In particular, the conditions on ∆x that rule out the existence of solutions of f (x) = y are sought. Relation (9) brings to

f C(y) + ∆x -y ∈ P C(y) + ∆x -y + I, (18)
which component-wise reads

f i C(y) + ∆x -y i ∈ P i C(y) + ∆x -y i + I i . ( 19
)
According to the mean value theorem (see Appendix A)

P i C(y) + ∆x = P i C(y) + ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x, (20)
for some suitable ξ i (y, ∆x) ∈ [0, 1]. Thus,

f i C(y) + ∆x -y i ∈ P i C(y) + ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x -y i + I i . ( 21
)
Using ( 12) yields

f i C(y) + ∆x -y i ∈ ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x + I i + Ji . (22)
i.e.,

f i C(y) + ∆x -y i ≥ ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x + inf(I i + Ji ) ≥ ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x -|I i + Ji |. ( 23
f i C(y) + ∆x -y i ≤ ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x + sup(I i + Ji ) ≤ ∇ x P i C(y) + ξ i (y, ∆x) ∆x • ∆x + |I i + Ji |. (24)
Select now the component ν that satisfies

|∆x ν | ≥ |∆x j | , ∀j. (25)
Hypotheses ( 13) and ( 14), and the arguments presented in Appendix B result in

f ν C(y + ∆x) -y ν ≥ d -k (n -1) ∆x ν -|I ν + Jν |, for ∆x ν ≥ 0 (26)
and

f ν C(y) + ∆x -y ν ≤ d -k (n -1) ∆x ν + |I ν + Jν |, for ∆x ν ≤ 0. ( 27
)
Based on (26), f ν C(y) + ∆x -y ν > 0 for any ∆x satisfying

∆x ν > |I ν + Jν | d -k (n -1) . (28)
Similarly, if

∆x ν < - |I ν + Jν | d -k (n -1) , (29)
then, from (27), f ν (C(y) + ∆x) -y ν < 0.

Consequently, there cannot be solutions of the problem f (x) = y, x ∈ D, for ∆x outside of the interval J = [-a, a] n , where a is defined as

a = max i=1,...,n |I i + Ji | d -k (n -1) , (30)
and all possible solutions must lie in the Taylor model C(y) + J .

Problem Formulation

Impulsive interplanetary transfers are usually designed in the frame of the patched-conics approximation. The hypotheses behind this approach can be listed as

(1) The trajectory is split into a sequence of two-body problems for which analytical solutions are available (conic sections) (2) The planets' spheres of influence are infinitely large in planetary reference frame but points in the heliocentric frame (3) Hyperbolic passages occur instantaneously (4) Different conic sections are patched together with instantaneous velocity changes delivered by the propulsion system. Engineering Optimization 11

For each arc, the dynamical model is described by Kepler's gravitational law

ṙ = v v = - µ r 3 r (31)
in which µ is the gravitational parameter of the central body and r = (x, y, z) and v = (v x , v y , v z ) the spacecraft position and velocity vectors, respectively. The solution to ( 31) is analytical and it is given by the conic section

r = p 1 + e cos ϑ , (32)
in which p is the semilatus rectus, e the eccentricity, and ϑ the true anomaly [START_REF] Kaplan | Modern spacecraft dynamics and control[END_REF]. If e = 0, the conic is a circle; if 0 < e < 1, the conic is an ellipse; if e = 1, the conic is a parabola; and, if e > 1, it is a hyperbola.

The transfer problem of interest is to bring a spacecraft from a circular orbit around a starting planet, typically the Earth, to a circular orbit around a target body. As the initial and the final orbit lie inside the planets' spheres of influence, the first and the last phases of the transfer are hyperbolic legs. These two legs are then connected outside the planets' spheres of influence by an elliptic heliocentric trajectory. The algorithmic flow for the computation of the objective function for a planet-to-planet transfer is

(1) Compute the position and the velocity of the planets at two given epochs t 1 and t 2 , where t 1 is the departure epoch and t 2 the arrival one (2) Determine the elliptic arc connecting the two planets in the transfer time ∆t = t 2 -t 1 (3) Compute the relative velocity at the planets, the two hyperbolic arcs, and the associated starting and arrival ∆v.

A different problem is associated with each of the former steps: the first requires the evaluation of ephemerides functions, the second is the solution of a two point boundary value problem known as Lambert's problem, and the last being simple orbital mechanics algebra. All these steps are classical astrodynamical problems, but they require particular attention when considered at the validated approach. For this reason each of these problems is examined separately in the following subsections. Note that throughout this paper the subscript 1 will denote starting conditions, 2 the arrival conditions. Furthermore, an Earth-Mars transfer will be adopted as reference planet-to-planet transfer.

Ephemerides

In order to design the transfer it is necessary to know planets' positions and velocities at given epochs. The ephemerides function used in this work is based on the ephemerides DE405 of Caltech's Jet Propulsion Laboratory (JPL). These have been obtained from a least-square fitting of previously existing ephemerides to the available observation data, followed by a numerical integration of a suitable set of equations describing the motion of the Solar system. The numerical integrations were carried out using a variable stepsize, variable order Adams method. The result of the integration is stored in form of interpolatory data (Chebyshev polynomials, each block of them covers an interval of 32 days). The DE405 ephemerides are valid from Dec 9, 1599 to Feb 1, 2200 [START_REF] Alessi | Efficient Usage of Self Validated Integrators for Space Applications[END_REF]). The internal reference system is the so-called J2000 coordinates. A detailed The evaluation of the positions and the velocities of planets within TM frame does not allow the use of any external code. In order to avoid this problem, a Taylor interpolation in time of planets' orbital parameters (a, e, i, ω, Ω, M ) obtained through JPL ephemerides has been carried out. In particular, a third order interpolation is selected to limit the interpolation error to the order of a few thousand km for the position and a m/s for velocities over the time windows of interest, an accuracy compatible with the preliminary optimization problems at hand. Conversion from orbit elements to cartesian quantities is performed as follows:

x = r cos(ϑ + ω) cos Ω -sin(ϑ + ω) cos i sin Ω y = r cos(ϑ + ω) sin Ω + sin(ϑ + ω) cos i cos Ω z = r sin(ϑ + ω) sin i v x = v -sin(ϑ + ω -γ) cos Ω -cos(ϑ + ω -γ) cos i sin Ω v y = v -sin(ϑ + ω -γ) sin Ω + cos(ϑ + ω -γ) cos i cos Ω v z = v cos(ϑ + ω -γ) sin i (33)
where the velocity v is

v = µ r - µ a , (34)
the true anomaly is related to the eccentric anomaly by

tan E 2 = 1 -e 1 + e tan ϑ 2 , ( 35
)
the flight path angle is obtained from

tan γ = e sin ϑ 1 + e cos ϑ , (36)
and the eccentric anomaly is related to the mean anomaly by Kepler's equation,

M = E -e sin E. (37)
Thus, based on the interpolated orbital parameters, the cartesian position and velocity of planets are computed as a function of time with only the major remaining difficulty of solving Kepler's equation in a validated way as described in Section 4.

Having the ephemerides evaluated in the TM framework gives the validated expansion of planets positions and velocities with respect to departure and arrival epochs. Figures 4 and5 are relative to Taylor model evaluations of Earth ephemerides in the interval epoch [3570,3580] MJD2000. The ratio between the interval widths (IW) and TM constant part is shown (in logarithmic scale) for both cases as a function of the expansion order. It is apparent that the IW reduces almost exponentially up to order 7 and that the 5-th order expansion guarantees good accuracy for the ephemerides evaluation.

∞ 1 v ∞ 2 v E v 1 v M v 2 r 1 r 2

Lambert's Problem

In Lambert's problem, the initial position, final position, and desired time for the transfer between the two positions are given. Solving Lambert's problem defines the orbital elements of the desired transfer orbit, allowing the calculation of the velocities at the initial and final position. A diagram for an Earth-Mars transfer is shown in Figure 6.

Lambert's theorem states that the time required for the transfer depends only on the semi-major axis, the sum of the two radii r 1 + r 2 , and the distance between the initial and final positions, i.e., the chord length c = |r 2r 1 | [START_REF] Kaplan | Modern spacecraft dynamics and control[END_REF]. The time required for the transfer can be written as

t 2 -t 1 = a 3 µ 2kπ + (E 2 -e sin E 2 ) -(E 1 -e sin E 1 ) . ( 38
)
The problem now is to find the correct values of a, E 1 , E 2 , and e that give the desired transfer time. As Lambert stated, however, the transfer time depends only on the three quantities mentioned earlier. The two radii and the chord length are already known from the problem definition. The semi-major axis is the only unknown parameter, and it follows All other variables in the equation should be known. This can be done [START_REF] Prussing | Orbital Mechanics[END_REF], or the transfer time can be written as a function of some other parameters such as p or ∆E. In our approach, based on Battin's algorithm [START_REF] Battin | An Introduction to the Mathematics and Methods of Astrodynamics[END_REF], the nonlinear equation to be solved is

log A(x) -log(∆t) = 0, ( 40
)
in which

A(x) = g(x) 3/2 α(x) -sin α(x) -β(x) + sin β(x) . ( 41
)
The functions α(x) and β(x) are related to x via the relations

sin 2 1 2 α(x) = s 2g(x) sin 2 1 2 β(x) = s -c 2g(x) , ( 42
) with g(x) = s 2(1 -x 2 ) , (43)
and the semiperimeter

s = (r 1 + r 2 + c)/2. ( 44
)
The logarithm ( 40) is adopted to facilitate the solution for x by means of standard Newton's iteration. Note that the relation between a and x is simply given by

a = s 2(1 -x 2 ) . ( 45
)
Once the nominal solution is available, Lambert's problem can be solved in a TM environment by applying the algorithm described in Section 4. The n-th order Taylor model of v 1 = v 1 (t 1 , t 2 ) and v 2 = v 2 (t 1 , t 2 ) are obtained as final output.

Objective Function Evaluation

The previous sections showed how to compute the conic arc that connects two planets for a given departure and arrival epoch (or transfer time). The procedure to compute the objective function including the departure and the arrival phases inside planets' spheres of influence is described in the following. After Lambert's problem has been solved, the outgoing and incoming velocities at planetary sphere of influences, v ∞ 1 and v ∞ 2 , are known, as shown in Figure 6, by and

v ∞ 1 = v 1 -v E (
v ∞ 2 = v 2 -v M . ( 47
)
As mentioned earlier, the goal of the interplanetary transfer is to bring the spacecraft from a circular orbit around the Earth to a circular orbit around the target planet. Furthermore, an altitude of 200 km is assumed for both orbits. Thus, two hyperbolic arcs can be designed, whose pericenters lie on the initial and final circular orbits of radius r c i , with i = 1, 2 as illustrated in Figures 7 and8.

The energy density for the two hyperbolic arcs is

E i = 1 2 v ∞ i 2 , ( 48
)
and the velocity at the pericenter is computed by solving

1 2 v p i 2 - µ r p i = E i (49)
for v p i . The circular velocities associated to the circular orbits are given by

v c i = µ/r c i . ( 50
)
As a result, the objective function is computed by

∆v = ∆v 1 + ∆v 2 = (v p 1 -v c 1 ) + (v p 2 -v c 2 ). ( 51
)
As spacecraft must be accelerated at departure and slowed down at arrival, each of the former ∆v is positive definite. objective function is included in the TM representation and the IW properly scales with the expansion order. Note that minimizing ∆v is equivalent to minimizing the propellant mass as, within the hypothesis of impulsive maneuvers, the propellant mass fraction is given by the well known rocket equation

f = 1 -e -∆v Ispg 0 , (52)
in which I sp is the engine specific impulse (typical value for chemical engines is in the range 250-400 s) and g 0 is the gravitational acceleration of the Earth at sea level.

Test Cases

Earth-Mars (EM) and Earth-Venus (EV) transfers are considered as test cases. These problems are bidimensional and they involve two ephemerides evaluations and one solution of the implicit equation representing Lambert's problem for each evaluation of the objective function. The departure epoch t 1 and time of flight ∆t are chosen as optimization variables in order to deal with a rectangular search domain. Over a sufficiently large search domain, determining the global optimum is difficult as the objective function shows a roughly periodic behavior, several bands of near-minima, and among these a large number of comparable local minima, as clearly highlighted in Figure 11.

One fundamental problem of the optimization at hand is the presence of discontinuity lines of the objective function within the search domain. Such discontinuities correspond to the so-called transitions from the "short-way" to the "long-way" solutions of Lambert 's problem, and vice-versa (Bernelli et al. 2007). A geometrical overview of the problem is reported in Figure 12 for a sample transition from the short-way to the long-way solution.

Transfer trajectories for an Earth-Mars transfer are plotted for a fixed time of flight of 320 days, moving from the left side of the discontinuity to the right side. On the left side of the discontinuity the short-way solutions are selected by the Lambert solver. Moving toward the right side, the orbital plane inclination of the transfer trajectories tends to increase. The discontinuity occurs when the transfer trajectory is exactly perpendicular to the ecliptic. Just after the occurrence of the discontinuity, in order to keep dealing with prograde solutions of the Lambert problem, the long-way solution is suddenly selected.

Corresponding to the previous transition, a plot of the z component of v ∞ 1 and v ∞ 2 is reported in Figure 13. The increase of the absolute values of v ∞ 1,z and v ∞ 2,z around t 1 = 1317 MJD2000 is associated to the increasing difference between the inclinations of the planetary orbital planes and the transfer trajectory; whereas the discontinuity is due to the switch from the short-way to the long-way solution. The discontinuity then affects all the algebraic operations ( 48)-( 51) required for the evaluation of the objective function.

In terms of the coding of the Lambert solver and the subsequent evaluation of the objective function, the discontinuity represents a transition from one branch of an if statement to another such branch. Unfortunately, this is a situation that cannot be handled in a direct fashion by self-verified software, since it would require the separate treatment of all solutions that fall into one branch of the if statement, and all other solutions that fall into the other branch.

In the discussion further below, some automatic methods that deal with this problem are outlined. However, in practical situations the problem is of reduced significance, since in a neighborhood of the transition, large values of ∆v necessarily appear, which automatically lead to a significant increase in the objective function that immediately precludes such regions. In fact, it is rather straightforward and common practice to perform simple estimates that determine an exclusion region around the discontinuities which can be rigorously excluded a priori without further work [START_REF] Izzo | Search space pruning and global optimisation of multiple gravity assist spacecraft trajectories[END_REF][START_REF] Bernelli-Zazzera | Global Trajectory Optimization: Can we Prune the Solution Space when Considering Deep Space Maneuvers[END_REF]. 46) and ( 47)) at the transition from the short-way to the long-way solution. [323.38, 323.48] [133.37, 133.47] 6.1. A Priori Pruning of Search Space Near Discontinuity

In the first approach to the problem, this practice is followed and common estimates are used to eliminate from the search space regions that cannot contain the optimum.

In order to conveniently characterize the regions of search space that cannot contain a minimum, the space is first discretized with a grid of equidistant boxes, and then all those boxes for which a simple a priori estimate shows that a minimum of the search space cannot occur are eliminated. For simplicity, subdomains are chosen that line up with the coordinates; but with limited extra effort, it is also possible to choose other objects, such as a discretization into parallelograms that more naturally follow the line of discontinuity. As long as this pruning is done in a rigorous way, which is not difficult because of the inherently large values of the objective function near the discontinuity, this leads to a validated strategy that does not impose significant burden, but eliminates the branching problem in the evaluation of the objective function, since now each subdomain of interest can be safely decided to fall into only one side of the branch.

The Earth-Mars case is now studied in the search space [3000, 4000] × [200, 400] MJD2000 × days subject to a priori pruning as shown in of Figure 14(a). Specifically, there is a discretization of boxes of 20 days of side, and boxes whose midpoint transfer radii are almost aligned (i.e., arcsin r1×r2 r1r2 < 10 deg) are pruned away. COSY-GO is then fed with 447 remaining boxes and run adopting a 5-th order expansion and setting the minimum box side to 0.1 day. [150, 600] [80, 350] When the larger search space defined in Table 2 is considered together with the same domain discretization and pruning criteria, the initial number of feasible boxes rises to 4899. In this case the number of steps required to compute the same global optimum enclosure is 25931 with a computational time of 715.2 s.

It is illuminating to observe the details of the optimizer's box splitting, reduction and elimination sequence. For this end, Figure 14 shows some zoom pictures in the optimizer's sequence of box splitting, reduction and elimination. Figure 14

Taylor & Francis and I.T. Consultant reduction tools; for example, the box covering the middle of the picture gets reduced by about a factor of five horizontally and a little more than two vertically, before being eliminated. The volume reduction of the box to the left of the middle box is even more dramatic; while a box to its right was first reduced only by small amount by LDB, and then subsequently split into two equal parts, before being eliminated.

The behavior near the global minimum is also rather revealing. Figure 14(c) shows the neighborhood of the minimizer at a scale of about 16 days by 14 days. It is seen that boxes at the bottom are reduced in domain by LDB before being eliminated, and boxes to the left and right are being significantly reduced. The last picture of the sequence shows the situation near the minimizer at a scale of about 14 by 7 days, again showing the interplay of domain reduction and box elimination.

Overall, a quantitative analysis of the active phase space region and the elimination pattern shows that boxes are eliminated with an average size of about 100 day squared, which must be considered rather favorable considering the fine detail in the objective function on the one hand, and its complexity involving a verified solution of an implicit equation on the other. It is also worthwhile to mention that there is no cluster effect problem visible, i.e. boxes near the minimizer follow a rather favorable pattern.

The search space of Table 2 is considered for the transfer to Venus. The same pruning technique adopted for the transfer to Mars is used, but on a discretization of boxes of 10 days. The optimizer is then fed with 12056 boxes. The number of iterations required is 88556 with a computational time of 2399.1 s. The enclosure of the solution is given in Table 1.

It is worth mentioning that if the departure and arrival legs are neglected (i.e., the objective function is simply computed by ∆v

= v ∞ 1 + v ∞
2 ) the enclosures of the minimum ∆v for the transfers to Mars and Venus are [5.659527, 5.659531] and [6.062911, 6.062919] km/s, respectively. These values compare well with the 5.6695 and 6.0638 km/s reported by [START_REF] Liu | The Application of Evolution-Branching Algorithm on Earth-Mars Transfer Trajectory[END_REF] and [START_REF] Bernelli-Zazzera | Global Trajectory Optimization: Can we Prune the Solution Space when Considering Deep Space Maneuvers[END_REF], which were obtained with different models and optimization methods.

Automated Pruning of the Search Space Near Discontinuity

It is worthwhile to study in detail the possibility of automating the pruning of the search space near the discontinuity lines. A first solution that is on the one hand somewhat simplistic, but on the other hand nevertheless allows for a detailed search of the parameter space, is to simply utilize the failure mode of the Taylor model implementation of COSY Infinity. Specifically, if there is a situation encountered where a Taylor model cannot be evaluated, an error flag is set instead, and all subsequent operations involving the erroneous Taylor model are themselves set to be erroneous. The COSY-GO verified global optimizer reacts to such a flag by simply splitting the relevant box of interest and trying again.

The advantage of this approach is its simplicity and rigor. It is certainly fully automatic, and any computation involving the (many) occurring erroneous Taylor models takes only minimal computational effort for propagating the corresponding error flag. On the other hand, the method leads to a simplistic box halving scheme that has the additional undesirable property that usually at least one of the split boxes is retained, namely the one containing the discontinuity line in the first place. This leads to a substantial increase in the number of boxes that need to be studied, in particular if the discontinuity regions are of extended size, as they are in the problem at hand. This approach is applied to optimize both Mars and Venus transfers on the search do- Engineering Optimization 21

Table 3. EM and EV transfers optimal solution enclosures with the automated pruning approach.

EM EV ∆v km/s [5.6974154, 5.6974160] [7.079469, 7.079472] t 1 MJD2000 [3574.28, 3574.38 ] [2715MJD2000 [3574.28, 3574.38 ] [ .69, 2715.77] ∆t day [323.41, 323.48] [133.33, 133.46] main of Table 2. Also for these optimization runs, a 5-th order expansion and a minimum box side of 0.1 day are chosen. The enclosure of the optimal ∆v and the corresponding optimal time variables are given in Table 3. In practice, the price to pay for obtaining the validated optimal solution in this simplistic manner is the computational time: the optimization processes require 202789 and 476425 iterations, corresponding to computational times of 7474.8 and 17557.1 s (on the same machine used before). Thus, this simplistic treatment of the pruning problem of the minimization leads to an increase of computational time by an order of magnitude.

In authors' opinion, the best method to attack the problem at hand in a fully automated fashion is the following. In case a subdomain is detected to contain a discontinuity, it is attempted to construct a sufficiently thin Taylor model "band" that contains the discontinuity. This can be readily achieved by utilizing an implicit function solver such as the one discussed in Section 4, but instead of using it to solve for the solution of Lambert's problem, use it to solve for the condition that produces the discontinuity. The resulting Taylor model including the discontinuity could then be broadened so that it can be safely eliminated.

Conclusion

The paper showed a first application of the Taylor-model-based global optimizer COSY-GO to space trajectory design. Simple impulsive planet-to-planet transfers were used as test cases, showing the efficacy of the method in finding the enclosure of the global minimum of the objective function. It is important to stress that the proposed optimization technique is deterministic, as opposed to all the available stochastic global optimizers, and it delivers the rigorous enclosure of the problem's global optimum. When used with a suitable a priori pruning technique to avoid undetermined regions of the objective function, favorable computation times can be achieved. On the other hand, fully automated detection of the undetermined regions requires further work; while it is possible to do this in an automated way based on arithmetic error traps, this approach leads to an increase of computational expense of one order of magnitude. Further work will focus on a more efficient representation of discontinuity bands arising from the Lambert solving process and the determination of sufficiently wide bands around them that the remaining active region always contain unique and continuous solutions. Future effort will be devoted to the reformulation of the problem in order to minimize the troubles associated to the discontinuities. This step is considered mandatory for the extension of the method to the rigorous global optimization of impulsive transfers that include multiple gravity assists and deep space maneuvers.

F

  Third order enclosure of the solutions of \mbox{$\sin (x) = y$} over $ D $.

Figure 1 .

 1 Figure 1. Third order enclosure of the solutions of sin(x) = y over D.

Figure 2 .Figure 3 .

 23 Figure 2. Accuracy of the enclosure of the solutions of sin(x) = y over D: comparison between different orders.

  Taylor & Francis and I.T. Consultant description about how these ephemerides are obtained is given by[START_REF] Standish | Detailed description of the JPL planetary and lunar ephemerides[END_REF].

Figure 4 .

 4 Figure 4. Taylor Model relative IW as a function of expansion order for Earth position.

Figure 5 .

 5 Figure 5. Taylor Model relative IW as a function of expansion order for Earth velocity.

  v

Figure 6 .

 6 Figure 6. Lambert's problem for Earth-Mars trajectories and starting and arrival velocity triangles definition.

  and I.T. Consultant that it is possible to write the transfer time simply as a function of the semi-major axis ∆t = f (a).(39)

Figure 8 .

 8 Figure 7. Sketch of Earth departure phase.

Figure 9 Figure 9 .

 99 Figure 9. Taylor Model evaluation of the objective function for a Earth-Mars transfer.

Figure 10 .

 10 Figure 10. Objective function Taylor Model relative IW as a function of the expansion order.

Figure 11 .

 11 Figure 11. Objective function landscape for Earth-Mars transfer.

Figure 12 .Figure 13 .

 1213 Figure12. Geometrical overview of the transition from the short-way to the long-way solution.

  (d) Zoom on box splitting, reduction and elimination sequence close to the global minimum.

Figure 14 .

 14 Figure 14. Analysis of Earth-Mars transfer domain reduction sequence.

  (b) shows boxes just at the point of being eliminated, at a location away from any local minimum. It is worth noting how boxes get shrunk significantly by the action of the LDB and QFB domain

Page 27 of 42 URL

	February 11, 2011 18:4 Engineering Optimization PaperGO-rev4
	Engineering Optimization
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

: http:/mc.manuscriptcentral.com/geno Email: A.B.Templeman@liverpool.ac.uk

Table 1 .

 1 EM and EV transfers optimal solution enclosures with the a priori pruning approach.

			EM	EV
	∆v km/s	[5.6974155, 5.6974160] [7.079470, 7.079472]
	t 1	MJD2000 [3574.30, 3574.37]	[2715.69, 2715.78]
	∆t day

  The optimization process requires 2233 steps and only 64.3 s of CPU time on a Xeon Dual Core 2.66 GHz laptop platform and ends with three remaining boxes, whose enclosure is given in Table1.

	February 11, 2011 18:4 Engineering Optimization PaperGO-rev4
			Engineering Optimization	19
		450			345
		400
					340
		350
	∆ t [AU]	300		∆ t [day]	335
		250
					330
		200
		2800 3000 3200 3400 3600 3800 4000 4200 150 t 1 [MJD2000] F (a) Earth-Mars initial list of boxes. o r 328	3570 (b) Box splitting, reduction and elimination se-3575 3580 3585 3590 3595 325 3600 t 1 [MJD2000] quence by LDB and QFB. 328
	∆ t [day]	324 326 322 314 316 318 320	P e r e	∆ t [day]	326 327 322 323 324 325
	quence close to the global minimum. (c) Box splitting, reduction and elimination se-e R 3566 3568 3570 3572 3574 3576 3578 3580 3582 312 t 1 [MJD2000] 321	3568 3570 3572 3574 3576 3578 3580 3582 t 1 [MJD2000]
				v i e
						w
						O n l
						y

Table 2 .

 2 EM and EV transfers search space definition.

		EM	EV
	t 1	MJD2000 [1000, 6000] [1000, 6000]
	∆t day

  ∈ D, x + ∆x ∈ D, and some suitable ξ i (x, ∆x) ∈ [0, 1]. For the sake of clarity, define h i = x + ξ i (x, ∆x) ∆x ∈ D, and rewrite (59) asg i (x + ∆x) = g i (x) + ∂g i ∂x i (h i ) ∆x i +Suppose now that the components of ∇ x g i can be bounded as

	February 11, 2011 18:4 Engineering Optimization PaperGO-rev4 February 11, 2011 18:4 Engineering Optimization PaperGO-rev4 February 11, 2011 18:4 Engineering Optimization PaperGO-rev4 February 11, 2011 18:4 Engineering Optimization PaperGO-rev4
	for x j =i	∂g i ∂x j	(h i ) ∆x j .	(60)
	o r F o r for all x ∈ D, with F o r F o r	∂g i ∂x i ∂g i (x) ≥ d > 0, ∀i ∂x j (x) ≤ k, ∀j = i,			(61)
	P e P e P P Select the component ν that satisfies d -(n -1) k > 0. e e			(62)
	Thus, from (60),	e r e r |∆x ν | ≥ |∆x j | , ∀j. e r e r			(63)
	R e R e R R g (64) e e
		v i e v i e v i e v i e
		w w w w
			O n l O n l O n l O n l
				y y y y

ν (x + ∆x) = g ν (x) + ∂g ν ∂x ν (h ν ) ∆x ν + j =ν ∂g ν ∂x j (h ν ) ∆x j ≥ g ν (x) + d ∆x ν -(n -1) k ∆x ν , for ∆x ν ≥ 0 ≤ g ν (x) + d ∆x ν -(n -1) k ∆x ν , for ∆x ν ≤ 0.

URL: http:/mc.manuscriptcentral.com/geno Email: A.B.Templeman@liverpool.ac.uk Engineering Optimization

Taylor & Francis and I.T. Consultant

Acknowledgments

The authors are grateful to Alex Wittig for the helpful discussions on the validated solution of implicit equations and to Ravi Jagasia for various fruitful discussions as well as a careful reading of the manuscript.

Appendix A

Given a scalar function g(x) : R n → R, with g continuous and differentiable, the aim of this Appendix is to show that

for some suitable ξ(x, ∆x) ∈ [0, 1].

Consider the continuous and differentiable function G(t) : [0, 1] → R defined as

The application of the mean value theorem for scalar functions to (54) yields

for some ξ(x)

Altogether, ( 55), (56), and (57) give

Appendix B

Given a vector-valued function g(x) : D ⊂ R n → R m , with g i continuous and differentiable on D for all i ∈ 1, . . . , m. According to the mean value theorem for multivariate functions reported in Appendix A