Abstract
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.
1. 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 problem 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 near-optimal 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.
2. Colliding Bodies Optimization Algorithm
Colliding Bodies Optimization is a metaheuristic algorithm developed by Kaveh and Mahdavi [11–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.
2.1. Basic CBO Formulation
Each search agent is modelled as a body with mass and velocity. The initial position of the th body is randomly provided in a -dimensional search space set by the user: where 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 (also called particle or object), its mass is defined as follows: where is the cost function value of the th particle and , which must be an even number, is the total number of bodies used in the optimization process (the population size). The 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:
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: 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: where 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.
| ||||||||||||||||
2.2. 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’ positions, the ECBO executes the following crossover instruction: where is the th variable of the th CB, randomly selected; are the lower and the upper bounds of the th variable; PRO is the crossover probability that must be set by the user between [0,1]; is a random number uniformly distributed within [0,1], automatically generated for each particle, as well as . If 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.
| ||||||||||||||||||
3. 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(3)The spacecraft’s initial position and velocity are the same as the Earth’s(4)The thrust magnitude of the spacecraft is constant
The spacecraft’s equations of motion are written in polar coordinates: where is the position vector; and are, respectively, the radial and the horizontal velocity of the spacecraft, is the Sun’s gravitational parameter, is the thrust magnitude; is the spacecraft mass, and is the thrust pointing angle relative to the local horizontal (Figure 1).
The control vector is . The thrust-to-mass ratio is as follows [7]: where is the initial spacecraft mass, is the exhaust velocity, and is the thrust-to-mass ratio at the initial time . A normalized set of units is as follows: (i)Distance unit is the mean radius of the Earth’s orbit(ii)Time unit , such that
The objective function to minimize is . 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 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: . 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 (), the same proportionality holds at any time. By writing the control as a function of the proportional set, it is possible to demonstrate that it coincides with the 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 . 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 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.
3.2. 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 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(2)The spacecraft’s initial position and velocity are the same as the Earth’s(3)The only attracting body is the Sun
A polar inertial reference frame is used, and the equations of motion for the solar sail are as follows: where 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. and 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); is the characteristic acceleration of the sail, and the terms , , and 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 , , and is considered [16].
The minimum time trajectory is obtained with an indirect approach. The Hamiltonian function is as follows: where are the time-dependent costate variables. The Euler-Lagrange equations are as follows:
The optimal control will be obtained in the form , but it is not possible to find an explicit solution in this form. The optimal steering law can be approximated by [17]: where
The equations of motion have the following boundary conditions:
There are four unknown parameters: , 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 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.
4. 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 is the position vector of the spacecraft and is the thrust vector. The mass of the spacecraft is indicated as and its time varying law is , where is the initial spacecraft mass and is the mass flow ratio. This set of equations is conveniently normalized using this set of units: (i)Distance unit is the mean radius of the Earth’s orbit(ii)Time unit , such that (iii)Mass unit
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 with variable duration , with . During each arc, the angles and represent the thrust vector direction in the local reference frame 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 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 , we need to minimize the following arrival errors as well: where subscript indicates the spacecraft state and 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 and (that weight of the different dimensionless contributes in the cost function) are 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 cost function evaluations . Each simulation will stop if 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 . A success occurs when the accuracy condition of the stopping criterion is met. Other performance indices are minimum and average values of transfer time and number of function evaluations . In addition, the dispersion around the mean value of transfer time will be shown.
4.1. Rendezvous with Mars
Here, the outer planet Mars has to be encountered by a spacecraft with the following characteristics: (i)Thrust (ii)Specific impulse (iii)Initial mass
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: (1)(2)(3)(4)
The number of arcs chosen is ; therefore, there are 13 variables. The coefficients of the cost function are so established: . The maximum number of function evaluations selected for this problem is 20000. It is used a basic version of the CBO by employing 50 particles to investigate the search space.
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.
4.2. Rendezvous with Venus
Here, the inner planet Venus must be encountered by a spacecraft with the same characteristics as the previous one.
The search space, in terms of angles, arc duration, and departure time, is defined as follows: (1)(2)(3)(4)
The number of arcs chosen is ; therefore, there are 16 variables. After some preliminary tests, the coefficients of the cost function are so established: . 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.
4.3. 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: (i)Thrust (ii)Specific impulse (iii)Initial mass
The search space is as follows: (1)(2)(3)(4)
The number of arcs chosen is resulting in 34 variables. The coefficients are . 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 is very 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 . 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.
5. 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 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.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.