Optimization of Interplanetary Trajectories Using the Colliding Bodies Optimization Algorithm

In this paper, a recent physics-based metaheuristic algorithm, the Colliding Bodies Optimization (CBO), already employed to solve problems in civil and mechanical engineering, is proposed for the optimization of interplanetary trajectories by using both indirect and direct approaches. The CBO has an extremely simple formulation and does not depend on any initial conditions. To test the performances of the algorithm, missions with remarkably different orbital transfer energies are considered: from the simple planar case, as the Earth-Mars orbital transfer, to more energetic ones, like a rendezvous with the asteroid Pallas.


Introduction
An important aspect of a space mission design is the obtaining of the nominal optimal trajectory, traditionally the one with minimum transfer time or maximum payload mass. Trajectory optimization is an old problem, its origins date back to the ancient Greeks, but its rigorous mathematical formulation, as an optimal control problem, arrives only with Pontryagin in the mid-1900s [1]. Optimal control is an issue concerning the determination of the inputs into a dynamical system that optimize (i.e., minimize or maximize) a specified performance index while satisfying several constraints [2]. These constraints can be differential, as the equations of motion, or algebraic, as departure, mid-course, and arrival constraints. Because of the complexity of most applications, optimal control problems are chiefly solved numerically. Numerical methods adopted are divided into two major classes: indirect methods and direct methods. In the former, the original problem is transcribed into a multiple-point boundary-value problem that is solved to determine candidate optimal trajectories, and the optimization phase consists in finding the optimal set of costate variables. In a direct method, the state and/or control of the optimal control prob-lem is discretized, and the problem is transcribed into a nonlinear optimization problem [2]. Both approaches lead to a parametric optimization where a set of optimal parameters must be found. This optimization has often been conducted by means of gradient-based research (e.g., Newton-Raphson-based algorithms), but in the last decades, a new kind of optimization procedures has been proposed and developed: metaheuristics [3,4]. The goal of metaheuristics is to efficiently explore the search space looking for nearoptimal solutions. The fundamental characteristics are that they are problem independent and they possess a nondeterministic nature, useful to escape from local optima. This is achieved by either allowing worsening moves or generating new starting solutions for the local search in a more "intelligent" way than just providing random initial solutions. This stochastic nature is not employed blindly, but in an intelligent, biased manner, and is what truly differentiates them from gradient-based techniques, which are deterministic and strongly dependent on an initial guess of the solution [5,6]. Gradient-based algorithms are largely employed in all fields of engineering, including space trajectories optimization. However, in recent years, metaheuristic algorithms have been increasingly adopted, especially in preliminary analysis of trajectories. Two of the most used algorithms, in this domain, are Particle Swarm Optimization [7,8] and Differential Evolution [9,10].
This paper is focused on the optimization of space trajectories using the Colliding Bodies Optimization algorithm. This is a novel population-based metaheuristic inspired by the one-dimensional collision theory between bodies, where each candidate solution being considered as a body with mass. CBO utilizes a simple formulation to find extremals of functions and does not depend on any internal parameter [11]. In Section 2, CBO and its enhanced version will be described briefly. Test cases using indirect methods are studied in Section 3, while those using direct methods are studied in Section 4. Conclusions will be given in Section 5.

Colliding Bodies Optimization Algorithm
Colliding Bodies Optimization is a metaheuristic algorithm developed by Kaveh and Mahdavi [11][12][13], inspired by the one-dimensional collision theory. There are two versions of this algorithm, a basic one [11] and an enhanced one [12], that improves the basic version by means of a sort of Elitism and Crossover. In the next two sections, some basic statements are reported while details can be found in the cited references.

Basic CBO Formulation.
Each search agent is modelled as a body with mass and velocity. The initial position of the ith body is randomly provided in a j-dimensional search space set by the user: where rand is a random number between 0 and 1. A collision occurs between two bodies, and their positions, after the impact, are updated based on the one-dimensional collision laws [11,13]. Given the body X k (also called particle or object), its mass is defined as follows: where J k is the cost function value of the kth particle and n, which must be an even number, is the total number of bodies used in the optimization process (the population size). The n colliding bodies (CBs) are sorted into ascending order, according to their objective function values, and then divided into two equal groups: Stationary Objects (the lower half) and Moving Objects (the upper half). Objects of the MO group collide against members of the SO group to improve their position and push stationary objects towards better positions. In particular, the colliding pairs are established according to the ascending order with respect to the objective function. Hence, for instance, the best moving particle collides with the best stationary one. Bodies' velocities before the collision are assigned as follows: Moving bodies As many other metaheuristic algorithms, velocities are not defined as the derivative of the position with respect to time, but they are expressed as displacements in the search space. According to the colliding bodies' theory, velocities after the collision are calculated as follows: Moving bodies where ε is the Coefficient of Restitution, defined as the ratio of the relative velocity between two bodies after and before the collision: This coefficient is assumed varying linearly between 1 and 0 during the optimization process, in order to ensure the balance between exploration and exploitation. After the calculation of the displacement, it is possible to determine new positions of the stationary and moving bodies as follows: Moving bodies : where rand is a uniformly distributed random vector in the range [-1,1]. This iterative scheme, performed on all the particles at each iteration, is repeated until a given stopping criterion is fulfilled. A Pseudocode 1 of CBO is reported.

Enhanced CBO Formulation.
The structure of the enhanced CBO (ECBO) algorithm is essentially the same as the basic CBO [12], with the difference that a Colliding Memory (CM) is introduced to save the best CBs' positions obtained so far. In fact, the positions stored in the Colliding Memory substitute the worst positions occupied by the current bodies. In this way, the best positions are remembered and there is no global worsening of the objective function from one iteration to another. The number of the best CBs' positions that are preserved, therefore the dimension of the Colliding Memory, is set by the user. Moreover, after the update of the CBs' 2 International Journal of Aerospace Engineering positions, the ECBO executes the following crossover instruction: if ran i < PRO, otherwise, where x ij is the jth variable of the ith CB, randomly selected; x j,min , x j,max are the lower and the upper bounds of the jth variable; PRO is the crossover probability that must be set by the user between [0,1]; ran i is a random number uniformly distributed within [0,1], automatically generated for each particle, as well as rand. If ran i is less than PRO, a crossover occurs. As PRO increases, the probability to perform a crossover increases. The Pseudocode 2 of the ECBO is reported.

Numerical Simulations: Indirect Methods
In order to analyze the performances of this optimization algorithm, a total of five study cases will be presented. Two of them are solved by using an indirect strategy, while the remaining three cases adopt a direct approach.
3.1. Optimal Earth to Mars Orbital Transfer. The problem is to reach the orbit of Mars departing from the Earth's orbit with the minimum transfer time by using a low thrust engine. This case has already been studied in literature with different techniques [7,14,15]. In this paper, the best trajectory is obtained by using the ECBO. As in the cited papers, the following hypotheses are established: (1) The orbits of the planets are coplanar and circular (2) The only attracting body is the Sun The spacecraft's equations of motion are written in polar coordinates: where r is the position vector; v r and v θ are, respectively, the radial and the horizontal velocity of the spacecraft, μ s is the Sun's gravitational parameter, T is the thrust magnitude; m is the spacecraft mass, and α is the thrust pointing angle relative to the local horizontal ( Figure 1). The control vector is uðtÞ = α. The thrust-to-mass ratio is as follows [7]: where m 0 is the initial spacecraft mass, c = 55:894 km/s is the exhaust velocity, and n 0 = 8:342 × 10 −4 m/s 2 is the thrust-tomass ratio at the initial time t 0 . A normalized set of units is as follows: (i) Distance unit DU = 149:5 ⋅ 10 6 km is the mean radius of the Earth's orbit (ii) Time unit TU = 5:018 ⋅ 10 6 s, such that μ s = DU 3 / TU 2 The objective function to minimize is J = t f . The desired terminal conditions are expressed by the vector ω as follows: To write the necessary conditions for optimality, the Hamiltonian function has to be defined: where the time-dependent set of costate variables ½p r , p v r , p v θ has been introduced. Furthermore, the costate differential equations are expressed as follows: From the Pontryagin principle, the optimal control law is as follows: The set of necessary conditions for optimality can be completed with the transversality condition referred to the final time: where λ is another set of adjoint variables concerning the terminal conditions; hence, it has the same dimensions as ω.
Rearranging Equation (17), the following condition is obtained: In order to find the optimal trajectory that satisfies the necessary conditions, the problem reduces to the determination of four parameters: ½p r ð0Þ, p v r ð0Þ, p v θ ð0Þ, t f . The CBO algorithm must find the solution exploring the following search space: It is necessary to introduce a cost function that links the optimization algorithm to the problem. This function always includes the quantity that must be minimized, and if there is any constraint to respect, the most popular approach is to include them in the cost function. Then, the objective function is defined here as follows: In the cost function, all the quantities (transfer duration and errors at final time) are dimensionless and the weights are chosen to scale them properly, in order to keep the balance between all terms throughout the simulation. The ECBO population is composed of 50 particles. The dimension of the Colliding Memory is set to 1 and crossover is not performed. The algorithm has to find the optimal set of costate initial values. The transversality condition can be neglected by the CBO, due to the homogeneity of the costate equations (Equation (15)). For this reason, if the CBO finds a set of initial costate variables that is proportional to the optimal one ( p  International Journal of Aerospace Engineering optimal control (Equation (16)). This means that the minimum time trajectory can be obtained also with an initial costate proportional to the optimal set. Nevertheless, the transversality condition will be violated, and by substituting the proportional set in Equation (18), it can be easily proved that Hðt f opt Þ = −b. In order not to increase the number of equality constraints in the cost function, the transversality condition is verified at the end of the optimization process, to guarantee the optimality of the trajectory. The optimal solution found by the CBO consists in a 192.6 days transfer ( Figure 2). Figure 3 shows the optimal thrust pointing angle αðtÞ and the costate evolution during the transfer trajectory. These results are in accordance with those in literature by using PSO [7] and are obtained quickly and easily thanks to the simplicity of the algorithm.

Optimal Earth to Mercury Orbital Transfer.
A minimum time transfer between the orbits of the Earth and Mercury is studied by means of a solar sail. The problem consists in determining the optimal steering law αðtÞ that minimizes the time of flight to reach Mercury's orbit [16].
As in the cited papers, a series of simplifying assumptions are made: (1) The relative orbital inclination of Mercury and Earth is neglected, and the orbits are considered circular 5 International Journal of Aerospace Engineering where r is the position vector and θ is the polar angle measured anticlockwise from the axis that connects the Sun to the Earth at the initial instant. v r and v θ are, respectively, the radial and the tangential velocity. α is the angle between the direction that connects the Sun to the solar sail and the thrust direction ( Figure 4); a c is the characteristic acceleration of the sail, and the terms b 1 , b 2 , and b 3 are the coefficients that represent the optical properties of the sail. A sail with an aluminium-coated front side and a chromium-coated back side with b 1 = 0:1728, b 2 = 1:6544, and b 3 = −0:0109 is considered [16].
The minimum time trajectory is obtained with an indirect approach. The Hamiltonian function is as follows: where ½p r , p θ , p v r , p v θ are the time-dependent costate variables. The Euler-Lagrange equations are as follows:   International Journal of Aerospace Engineering The optimal control will be obtained in the form α = αðp v r , p v θ Þ, but it is not possible to find an explicit solution in this form. The optimal steering law can be approximated by [17]: where α * p ≅ 2:5392 rad, α ≅ 0:008109α 6 p − 0:05474α 5 p + 0:1356α 4 p − 0:1266α 3 p + 0:08266α 2 p + 0:3038α p + 0:0008666 with α p ∈ 0, α * p h i : The equations of motion have the following boundary conditions: There are four unknown parameters: ½p r ð0Þ, p v r ð0Þ, p v θ ð0Þ, t f , and the optimal values are found in the following search space: To investigate the domain, a population of 50 particles is considered. As in the previous case, an ECBO that preserves the best position in the population at each iteration is employed. Therefore, the selected dimension of the Colliding Memory is 1 and the crossover is not considered. The cost function is written, as in the previous case:  Although the cost function has the same form as the previous one, coefficients must be chosen carefully and they are different for each problem. The characteristic acceleration is set to 0:25mm/s 2 and the resulting transfer lasts 2.8 years, in agreement with the results obtained in [16]. The optimal costate and α time history are shown in Figure 5. The CBO, also in this case, demonstrates a great ability in finding the optimal initial costate values within 6500 function evaluations with no parameter tuning at all.

Numerical Simulation: Direct Methods
In this case, the test cases proposed are minimum time rendezvous with three different celestial bodies: an outer planet, an inner planet, and a very inclined asteroid. The transfer is considered heliocentric and the mathematical model consists in a three-dimensional restricted two body problem, with the Sun as the attractor and the spacecraft with negligible mass. The equations of motion are as follows: where r = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðx 2 + y 2 + z 2 Þ p is the position vector of the spacecraft and T is the thrust vector. The mass of the spacecraft is indicated as m and its time varying law is m = m 0 − _ m ⋅ t, where m 0 is the initial spacecraft mass and _ m is the mass flow ratio. This set of equations is conveniently normalized using this set of units: A direct single shooting technique, with constant thrust magnitude, is proposed due to its simplicity and effectiveness [18]. The trajectory is divided into a finite number of arcs NA with variable duration Δt k , with k = 1, ⋯, NA. During each arc, the angles α k and ψ k represent the thrust vector direction in the local reference frame ðr, b θ,ĥÞ with respect to the inertial reference frame J2000 (Figure 6).
In a rendezvous problem, the departure time is another variable to compute during the optimization process. The number of arcs NA is selected after a preliminary study and chosen to reduce as much as possible the computational effort [19]. Globally, there will be three variables for each arc plus the departure time. The ephemerides of the planets are DE430 (provided by the NASA SPICE tool [20]). With the transfer time t f , we need to minimize the following arrival errors as well: where subscript s indicates the spacecraft state and t represents the target body. The desired final condition, hence, is that the arrival errors (Equation (30)) are zero. The objective function is defined as follows: The coefficients C t , C d and C v (that weight of the different dimensionless contributes in the cost function) are  Figure 6: Thrust vector control direction in the local reference frame and its orientation with respect to J2000. chosen after a preliminary study to properly balance the different terms, and their values will be presented for each test case. In the next section, three rendezvous missions will be discussed and the selected celestial bodies are Mars, Venus, and Pallas. The stopping criterion is composed of two conditions. The first one concerns the accuracy of the mission: each transfer will be considered successful, and hence, the relative simulation will stop, if the errors on position and velocity, described in Equation (30), go, respectively, below the mean radius of the target body and the threshold velocity of 20 m/s. The second term of the stopping criterion is a computational condition that imposes a maximum number of   9 International Journal of Aerospace Engineering cost function evaluations F eval . Each simulation will stop if F eval exceeds the threshold value set for the specific mission. If this threshold is reached, the algorithm automatically reinitializes itself with a different initial distribution of particles. Both versions of the CBO will be used. In order to perform a more trustworthy analysis, one thousand runs will be executed for each case. Each run departs from different initial distributions of particles in the search space because of the randomness of the initialization of the positions (Equation (1)). This is done to evaluate the average performances of the algorithm. The first performance index considered is the percentage of success ð% succ Þ. A success occurs when the accuracy condition of the stopping criterion is met. Other performance indices are minimum and average values of transfer time ðt f Þ and number of function evaluations ðF eval Þ. In addition, the dispersion around the mean value of transfer time σ t f will be shown.

Rendezvous with Mars.
Here, the outer planet Mars has to be encountered by a spacecraft with the following characteristics: (i) Thrust T = 300 mN (ii) Specific impulse I sp = 3000 s (iii) Initial mass m 0 = 1000 kg The search space, in terms of angles, arc duration, and departure time, must be chosen to let the optimization process start. The optimizer will look for the best solution out of the bounds here empirically defined: (3) Δt k ∈ ½10, 100 days   As can be seen from Table 1, the minimum time transfer, on average, lasts 302 days with a minimum of 288 days. The standard deviation is small; therefore, although the search space is wide, the CBO finds similar trajectories and this demonstrates the capability of the algorithm to explore effectively the search space. In Figure 7, the behavior of the cost function is plotted for all the successful runs (77.8%). It can be seen that the cost function decreases similarly for each run, which can be interpreted as a sign of robustness of the algorithm at changes in the initial conditions. The best trajectory found so far carries the spacecraft to Mars in 288 days with a final mass of 746.35 kg (Figure 8). Finally, the control law is plotted in Figure 9. The search space, in terms of angles, arc duration, and departure time, is defined as follows: (1) α k ∈ ½−180°, 0° (2) ψ k ∈ ½−90°, 90° (3) Δt k ∈ ½1, 90 days (4) t dep ∈ ½09/01/2019 − 09/01/2020 The number of arcs chosen is NA = 5; therefore, there are 16 variables. After some preliminary tests, the coefficients of the cost function are so established: C t = 0:01, C d = 100, C v = 50. In this problem, the number of maximum function evaluations is set to 50000. It is employed a basic CBO with 80 particles. The results are shown in Table 2. The ability of the CBO to both explore and exploit is confirmed here by the small dispersion of the duration transfer around the mean value.
The balance between exploration and exploitation phase can also be deduced by the almost linear slope of the cost function trend, plotted in Figure 10. The optimal control law and the optimal trajectory are shown in Figures 11 and 12.

Rendezvous with Pallas.
Pallas is an asteroid belonging to the asteroid belt that describes a very inclined orbit around the Sun (34.8°) with an eccentricity of 0.2305, making it a difficult body to reach. The probe chosen for this mission has the following characteristics: The search space is as follows: (1) α k ∈ ½0°, 180° (2) ψ k ∈ ½−90°, 90° (3) Δt k ∈ ½5, 250 days   The maximum number of function evaluations allowed for this problem is 100000. It is employed an ECBO with one particle in the Colliding Memory without considering the crossover. The number of particles is 120. Despite the biggest energetic difference between the Earth and target orbits, the strictest arrival tolerances, the least powerful thruster, and the highest number of variables, the CBO can find optimal solutions for this rendezvous. The success percentage of 23% (Table 3) can be addressed to the concepts just mentioned, yet the value of standard deviations of t min is very

13
International Journal of Aerospace Engineering small (3% of the mean transfer time), revealing that the discovered trajectories are very close to one another. This can be interpreted as a sign of the reliability of the algorithm. Figure 13, Figure 14, and Figure 15 represent, respectively, the optimal control law, the cost function behavior, and the optimal trajectory. 4.4. Discussion of the Results. It is worth mentioning that the direct approach, utilized in the last three cases, does not contemplate optimality conditions; therefore, there is not a definite optimal trajectory. Since there is no reference solution, a large number of simulations is needed to characterize the behavior of the algorithm. Due to the many degrees of freedom offered, the developed strategies allow many solutions to satisfy the arrival constraints. Although the range of possible transfer times is wide, the solutions found by the CBO are not equally distributed on all the possible values of t f . In Figures 16, 17, and 18, the transfer durations of the successful trajectories are grouped into histograms. For all the cases tested above, they are centered around a mean value, biased towards the minimum time transfer, with very small standard deviation values (equal or less than 3% of the mean value), reported in Tables 1, 2, and 3.

Conclusions
The study conducted in this paper reveals that the CBO is a valid algorithm, capable of solving a broad range of trajectory optimization problems. In the indirect cases, the CBO finds the optimal set of costate variables, discovering the optimal trajectory in a straightforward way. The direct cases were  14 International Journal of Aerospace Engineering made challenging on purpose, to test the effectiveness of the CBO thoroughly. Despite the increasing number of variables, the large domains, and the strict arrival tolerances, the algorithm succeeded in achieving satisfying trajectories for all the problems presented, exhibiting also a small dispersion around the suboptimal solution. Concerning the CBO performances, it is possible to state that, although it is straightforward in terms of computational effort and parameter settings, it shows great reliability and effectiveness in solving these kinds of problems. The remarkable advantage of CBO is that it is ready to use and it does not need an initial guess of the solution or insights of the problem. In addition, CBO has just one parameter to adjust (Elitism scheme); therefore, it does not need the fine-tuning preliminary operations required by metaheuristics in general.

Data Availability
The data used to support the findings of this study are included within the article.