Optimization of the Forcing Term for the Solution of Two-Point Boundary Value Problems

Abstract

We present a new numerical method for the computation of the forcing term of minimal norm such that a two-point boundary value problem admits a solution. The method relies on the following steps. The forcing term is written as a (truncated) Chebyshev series, whose coefficients are free parameters. A technique derived from automatic differentiation is used to solve the initial value problem, so that the final value is obtained as a series of polynomials whose coefficients depend explicitly on (the coefficients of) the forcing term. Then the minimization problem becomes purely algebraic and can be solved by standard methods of constrained optimization, for example, with Lagrange multipliers. We provide an application of this algorithm to the planar restricted three body problem in order to study the planning of low-thrust transfer orbits.

1. Introduction

We consider a two-point boundary value problem for the nonautonomous dynamical system in ℝ^𝑁: 𝑦^′(𝑡)𝑦(𝑡)(𝑡),𝑦()=𝑓+𝜑−1=𝑦₀,𝑦(1)=𝑦₁, (1) where 𝑓∈𝐶¹(ℝ^𝑁,ℝ^𝑁) and 𝜑∈𝐶([−1,1],ℝ^𝑁). We introduce a new numerical method for the computation of 𝜑 of minimal norm, among the functions representable with a truncated Chebyshev series, such that (1) admits a solution. In order to evaluate the performance of the algorithm, we apply it to a well-known problem, that is, the study of low-thrust orbits between the Earth and the Moon. We refer to [1] and references there for a discussion of the astrodynamical problem.

Our new approach is derived from the method introduced in [2, 3] to consider the dependence of parameters of the solution of hyperbolic equations. Such method heavily relies on a computer representation of algebras of functions, in such a way that the computer can use such functions in a straightforward and transparent way as if they were floating point numbers, reducing to a minimum the computations to be done explicitly by hand. The underlying ideas of this method stem from automatic differentiation algorithms, see [4–8] and recently have been used in the field of computer-assisted proofs, see [9, 10] for some examples, but its potentiality in optimization and control appears to have been overlooked. In Section 2 we describe the approach in a generic setting; in Section 3 we describe the specific-planar-restricted three body problem that we chose as an example; in Section 4 we present the results that we obtained with the three body problem; in Section 5 we discuss the features of the method and compare the Taylor and the Chebyshev representation.

2. The Generic Problem

Consider an initial value problem in ℝ^𝑛𝑦^′_𝑐(𝑡)=𝑓𝑡,𝑦(𝑡)+𝜑_𝑐(𝑡),𝑦_𝑐(−1)=𝑦₀, (2) where 𝑓 is an analytic function, 𝑐={𝑐_𝑖𝑗} is an 𝑛×𝑁 real matrix and the forcing term 𝜑_𝑐∶[−1,1]→ℝ^𝑛 depends on 𝑐 as follows: 𝜑_𝑐(𝑡)_𝑖=^𝑁−1_𝑗=0𝑐_𝑖𝑗𝑇_𝑗(𝑡), (3) where 𝑇_𝑗(𝑡) is the 𝑗th Chebyshev polynomial and (𝜑_𝑐(𝑡))_𝑖 denotes the 𝑖th component of 𝜑_𝑐(𝑡). Our first step is to compute explicitly the dependence of 𝑦_𝑐(𝑡) on 𝑐, for 𝑡∈(−1,1].

We choose a numeric method to solve the IVP; in order to focus on the issue of the parameter dependence, we choose the simplest method, for example, an explicit first-order Euler scheme with constant time step, and we consider the case 𝑛=𝑁=1. We show later that it is straightforward to adapt the method to other schemes, for example, a Runge-Kutta scheme with variable time step, and to the case 𝑛, 𝑁>1.

Choose 𝐾∈ℕ, let ℎ=2/𝐾 be the time step, set 𝑡_𝑘=𝑘ℎ−1 and 𝑌_𝑘(𝑐)=𝑦_𝑐(𝑡_𝑘) so that the Euler scheme reads 𝑌_𝑘+1(𝑐)=𝑌_𝑘(𝑐)+ℎ𝑓𝑡_𝑘,𝑌_𝑘(𝑐)+ℎ𝜑_𝑐𝑡_𝑘 (4) and 𝑌_𝐾(𝑐)=𝑦_𝑐(1). Now write 𝑌_𝑘(𝑐)=^+∞_𝑗=0𝑦_𝑘,𝑗𝑈_𝑗(𝑐), (5) where {𝑈_𝑗(𝑐)}_𝑗∈ℕ is a basis of some space of functions 𝑋 of our choice. In this paper we consider both the Taylor expansion for the set of analytic functions in a disc (i.e., 𝑈_𝑗(𝑐)=𝑐^𝑗) and the Chebyshev expansion for the set of continuous functions in [−1,1] (i.e., 𝑈_𝑗(𝑐)=cos(𝑗arccos𝑐)). Clearly, the method can be generalized to other expansions, for example, Fourier series or Legendre polynomials. Now we choose an order 𝑀 and we approximate 𝑌_𝑘(𝑐) with the truncated expansion ∑^𝑀_𝑗=0𝑦_𝑘,𝑗𝑈_𝑗(𝑐). Our aim is to compute the coefficients 𝑦_𝑘,𝑗 using an automatic algorithm. Since 𝑦₀(𝑐)=𝑦₀, then we have 𝑦_0,0=𝑦₀ and 𝑦_0,𝑗=0, when 𝑗≥1.

The core of the method consists in a generalization of the Taylor Models approach; for a detailed description of the technique with different expansions see [2, 3]. Here, we only recall the basic ideas. First, we observe that it is straightforward to implement on a computer an arithmetic of Taylor or Chebyshev polynomials of one variable. Using object-oriented programming and operator overloading, it is possible to define a class Fun which represents a Taylor or Chebyshev expansion and a set of methods which perform the basic operations and the computation of the elementary functions. More precisely, an object Fun is represented as the list of the coefficients. It is straightforward to implement a procedure that, given a scalar 𝛼 and two objects Fun  𝑇₁,𝑇₂, computes the object Fun corresponding to 𝛼(𝑇₁∗𝑇₂), where ∗ is either the addition or the multiplication. Then, given a polynomial 𝑝(𝑥), it is possible to compute the object Fun corresponding to 𝑝(𝑇₁), and since all analytic functions 𝑓(𝑥) can be approximated with polynomials, it is also possible to compute the object Fun that better approximates 𝑓(𝑇₁). A more algebraic-oriented way of expressing this idea is the following: let 𝑋 be the space of functions spanned by {𝑈_𝑗(𝑐)}_{𝑗∈{0,…,𝑀}} and let 𝑇∶𝑋→ℝ^𝑀+1 be the (invertible) map that extracts the coefficients. Our approach consists in lifting the basic arithmetics of 𝑋 to ℝ^𝑀+1. More precisely, we compute operators {⊕,⊗}∶ℝ^𝑀+1×ℝ^𝑀+1→ℝ^𝑀+1,⊙∶ℝ×ℝ^𝑀+1→ℝ^𝑀+1, (6) such that 𝑇(𝑓+𝑔)=𝑇(𝑓)⊕𝑇(𝑔), 𝑇(𝑓𝑔)=𝑇(𝑓)⊗𝑇(𝑔), and 𝑇(𝛼𝑔)=𝛼⊙𝑇(𝑔). Once this has been done, we can make extensive use of object-oriented programming and operator overloading in order to make this lifting completely transparent to the user.

The extension of this method to the case 𝑛>1 is straightforward. The case 𝑁>1 requires some additional work, since we have to consider a multivariable Taylor or Chebyshev expansion: 𝑌_𝑘(𝑐)=^+∞_{𝑗1,…,𝑗𝑁=0}𝑦_{𝑘,𝑗1,…,𝑗𝑁}𝑈_𝑗1𝑐₁⋯𝑈_𝑗𝑁𝑐_𝑁. (7) The sum of two of such expansions is done componentwise. The multiplication is less straightforward, but a very efficient algorithm for the case of the Taylor expansion has been introduced in [8], and it is not too hard to extend it to the Chebyshev (or some other) expansion. We refer to our programs for details.

With the class Fun consisting of the objects (7) and the functions for addition and multiplication in place, we can use a standard implementation of the algorithm (4), observing that, given 𝑌_𝑘, the methods of the class will take care of the computation of 𝑓(𝑡_𝑘,𝑌_𝑘), while 𝜑_𝑐(𝑡_𝑘) is explicitly given by (3). Extending this approach to better integration schemes, say Runge-Kutta or multistep, is straightforward: it suffices to implement the chosen algorithm using objects Fun and their methods. By implementing and running the integration scheme of choice with the class Fun, we obtain all the coefficients {𝑦_{𝑘,𝑗1,…,𝑗𝑁}} and therefore the polynomials 𝑌_𝑘(𝑐) for all 𝑘. Now consider the polynomial 𝑌_𝐾(𝑐): solving the boundary value problem (1) consists in computing 𝑐∈ℝ^𝑛𝑁 such that 𝑌_𝐾(𝑐)=𝑦₁. This amounts to solving 𝑛 algebraic equations in 𝑁 unknowns. If we wish ‖𝜑_𝑐‖ to be small (where ‖⋅‖ is some norm of our choice), we recall that when 𝑁>1 the system 𝑌_𝐾(𝑐)=𝑦₁ is underdetermined, therefore we can compute an explicit formula Φ(𝑐)=‖𝜑_𝑐‖ and exploit the additional degrees of freedom to set the following constrained minimization problem: find 𝑐∈ℝ^𝑛𝑁 such that Φ(𝑐)=min_{𝑌𝐾(𝑐)=𝑦1}Φ(𝑐). (8) This result can be achieved by standard constrained minimization algorithms, for example, Lagrange multipliers. Note that this method can be easily extended to solve other kinds of conditions at 𝑡=1: let 𝑚≤𝑛 and 𝑓∶ℝ^𝑛→ℝ^𝑚 some polynomial function. We can apply the same approach to find 𝑐∈ℝ^𝑛𝑁 such that Φ(𝑐)=min_{𝑓𝑌𝐾(𝑐)=0}Φ(𝑐). (9) To sum up, once the class Fun and some algorithm for solving initial value problems for ODE's has been implemented, boundary value problems with optimization of the forcing term are reduced to computing the constrained minumum of a polynomial, the constraint being represented by a polynomial equation or system of equations. We remark that, in principle, such polynomial constrained minimization problems can be solved exactly (up to floating point errors), but this does not imply that the boundary condition is solved exactly for the original problem, since we are truncating the expansion (7) and we are neglecting the truncation error. One needs to check a posteriori, for example, by solving (2) with standard floating point numbers and with the forcing term provided by (3), with 𝑐=𝑐 and computing the error 𝜀=|𝑦(1)−𝑦₁| or 𝜀=|𝑓(𝑦(1))|.

3. An Application: Low-Thrust Orbits in the Planar RTBP

In order to model the dynamics of a spaceship travelling from an orbit around the Earth to an orbit around the Moon, we consider the forced planar RTBP, that is the system: 𝑥̈+2𝑦̇=Ω_𝑥(𝑡)𝑥̇+𝑓𝑥̇²+𝑦̇²,𝑦̈−2𝑥̇=Ω_𝑦(𝑡)𝑦̇+𝑓𝑥̇²+𝑦̇², (10) where Ω=𝑥𝑥,𝑦²2+𝑦²2+1−𝜇𝑥+𝜇²+𝑦²+𝜇𝑥+𝜇−1²+𝑦²+𝐶. (11) We choose 𝜇≃0.0121506683, that is, the reduced mass of the Earth-Moon system, and 𝐶=−1.600172454916536, so that Ω(𝐿₁)=0, 𝐿₁ being the Lagrangian point between the primaries. The forcing term is given by a scalar function 𝑓(𝑡) to be determined times a unit vector directed as the velocity of the spaceship. In other words, we assume that the thrust is always parallel to the velocity. We show how to apply the method described in the previous section to the problem of optimizing the travel from an orbit 167 km above the Earth to an orbit 100 km above the Moon.

We choose a connecting orbit passing close to 𝐿₁. Note that the Jacobi function 𝐽𝑥,𝑦,𝑥̇,𝑦̇=𝑥̇²+𝑦̇²−2Ω𝑥,𝑦 (12) is approximately −60 on an orbit 167 km above the Earth and approximately −1.6 on an orbit 100 km above the Moon. Furthermore, in order to “cross” the region around the Lagrangian point 𝐿₁, it is necessary to raise it above 0. We build the connecting orbit as follows. We start from the Lyapunov orbit around 𝐿₁ with 𝐽=0.1 and approximate periodic orbits around the Moon and the Earth. Then we compute approximate connections from the Lyapunov orbit to the orbit around the Earth (backwards in time) and to the orbit around the Moon (forward in time). Choosing the Lyapunov orbit as a starting point is useful, if not necessary, because of its high instability: it is much harder to aim at the Lyapunov orbit than at a stable orbit around a primary. These orbits are obtained by applying a constant low thrust (i.e., 𝑓(𝑡)=𝑐, with small 𝑐), chosen by simple trial and error.

Once these rough orbits are computed, we apply the technique described above to compute an accurate correction so that the trajectory satisfies some required boundary conditions and at the same time the thrust is minimal. More precisely, we write (10) as a system of 4 first-order equations by setting 𝑣=𝑥̇ and 𝑤=𝑦̇, we write the modulus of the forcing term as 𝑓(𝑡)=𝑐+^𝑁−1_𝑗=0𝑐_𝑗𝑇_𝑗(𝑡), (13) we choose some initial value on the trajectory at the arbitrary time 𝑡=0 and we solve the equation with a fourth-order Runge-Kutta scheme up to some time 𝑡=𝑇, using objects Fun as scalars, the coefficients 𝑐_𝑗 being the independent variables of the objects Fun. Referring to Section 2, we have 𝑛=4 (the dimension of the system), we choose 𝑁=4 and 𝑀=10. We rescale the time in the interval [−1,1] and apply the method described in Section 2 to obtain (𝑥(𝑇),𝑦(𝑇),𝑣(𝑇),𝑤(𝑇)) as a function of {𝑐_𝑗}, expressed either as a Taylor series or a Chebyshev series.

Then we choose the condition to be satisfied at 𝑡=𝑇. We tested the following: find {𝑐_𝑗} such that (𝑥(𝑇),𝑦(𝑇),𝑣(𝑇),𝑤(𝑇)) is the same as in the case 𝑓(𝑡)=𝑐, but ‖𝑓‖ is minimal; find {𝑐_𝑗} such that the distance of the spaceship to the Earth or the Moon is assigned together with the modulus of the final velocity, that is, (𝑥(𝑇)−𝜇)²+𝑦(𝑇)−𝑅²=0 and (𝑣(𝑇))²+𝑤(𝑇)²−𝑉²=0 and ‖𝑓‖ is minimal; find {𝑐_𝑗} such that the distance of the spaceship to the Earth or the Moon is assigned together with the modulus of the final velocity, with the additional requirement that the final velocity is orthogonal to the line connecting the satellite to a primary. All these problem belong to the class of problem described in the previous section, and the polynomial constrained minimization problems have been solved by the Lagrange multipliers method.

4. Results

Figure 1 represents an example of an orbit computed by using the technique. Tables 1 and 2 collect some results obtained for a trajectory backward in time starting close to a Lyapunov orbit around 𝐿₁, more precisely at the point (𝑥(0),𝑦(0),𝑣(0),𝑤(0))=(0.73356888,−0.017228233,0.25064405,0.17220602) and ending close to an orbit around the Earth at time 𝑇=−1.7. The trajectory with a constant thrust 𝑐=0.2 ends at 𝑅=(𝑥(𝑇)−𝜇)²+𝑦(𝑇)=0.421852 and 𝑉=(𝑣(𝑇))²+𝑤(𝑇)²=0.876662.

Table 1 

Result enforcing 𝑅 and 𝑉.

Table 2 

Result enforcing 𝑅, 𝑉 and tangential trajectory.

Figure 1 

A complete orbit.

For both tables we computed the thrust of minimal norm necessary to end the trajectory with the values given in the columns 𝑅 and 𝑉. The results of Table 2 represent the solution with the additional constraint that the final velocity is orthogonal to the segment from the spaceship to the Earth, corresponding to an approximate circular orbit. We tested the accuracy of the computation as explained at the end of Section 2, that is, we repeated the computation of the orbit with standard floating point numbers, using the explicit expression of the forcing term, and we used the result of such computation to verify that the constraints were satisfied. The values 𝜀_𝑇 and 𝜀_𝐶 represent the mean square of the relative errors on 𝑅 and 𝑉, obtained with the Taylor and the Chebyshev computation, respectively, while 𝜀^⟂_𝑇 and 𝜀^⟂_𝐶 is |(𝑥(𝑇)−𝜇)𝑣(𝑇)+𝑦(𝑇)𝑤(𝑇)|, that is, the error on the orthogonality requirement. The columns ‖𝑓_𝑇‖ and ‖𝑓_𝐶‖ represent the norm of the forcing terms. It is clear from the tables that the method can achieve very accurate results, and it is also clear that the norm of the forcing term obtained with the Taylor and the Chebyshev computation is (almost) the same.

5. Conclusions

We discuss and compare the features of the two expansion that we tested: Taylor and Chebyshev. Generally speaking (see [2, 3]), the main advantages of the Taylor expansion consist in three features: the algorithms are simpler and faster, they provide directly the derivatives with respect to the parameters and furthermore, and, since one does not need to know the radius of convergence a priori, one can apply them and compute a posteriori the interval of validity of the computation. The Taylor expansion has also two main disadvantages: it does not provide uniform errors in the interval where the computation is performed and the radius of convergence may be bounded by poles in the complex plane. The Chebyshev expansion has opposite features: the error is uniform in the interval and the region of convergence is an ellipse as thin as necessary, therefore there are no problems caused by poles with nontrivial imaginary parts. On the other hand, the algorithm for the multiplication is much more complicated and much slower, it does not provide any direct computation of the derivatives and finally, since the Chebyshev polynomials are defined in [−1,1], one needs to choose a priori (via a suitable translation/rescaling) the interval where the parameter ranges, finding out only a posteriori if the approximation is acceptable. In the application considered in [2, 3], the Chebyshev expansion was a clear winner, due to the fact that, when considering the numerical approximation of shock waves and when trying to improve the resolution, the problem of poles becomes the main issue. Here, our judgment is the opposite: since the solutions of the ode's that we are studying are analytic functions, the issue of poles limiting the radius of convergence of Taylor series becomes negligible, while the disadvantages of the Chebyshev expansion, in particular the slower algorithms and the fact that the domain has to be chosen a priori, become very relevant. Additionally, we do not have evidence of a better performance in terms of accuracy. It is easy to find examples when either expansion performs better than the other one, but on average the accuracy is similar.